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ABSTRACT This paper reviews a class of generic dissipative dynamical 
systems called N-K models. In these models, the dynamics of N elements, 
defined as Boolean variables, develop step by step, clocked by a discrete 
time variable. Each of the N Boolean elements at a given time is given a 
value which depends upon K elements in the previous time step. 
We review the work of many authors on the behavior of the models, looking 
particularly at the structure and lengths of their cycles, the sizes of their 
basins of attraction, and the flow of information through the systems. In 
the limit of infinite N, there is a phase transition between a chaotic and an 
ordered phase, with a critical phase in between. 

We argue that the behavior of this system depends significantly on the 
topology of the network connections. If the elements are placed upon a 
lattice with dimension d, the system shows correlations related to the stan- 
dard percolation or directed percolation phase transition on such a lattice. 
On the other hand, a very different behavior is seen in the Kauffman net 
in which all spins are equally likely to be coupled to a given spin. In this 
situation, coupling loops are mostly suppressed, and the behavior of the 
system is much more like that of a mean field theory. 

We also describe possible applications of the models to, for example, genetic 
networks, cell differentiation, evolution, democracy in social systems and 
neural networks. 
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1 Introduction 

In this review, we describe the dynamics of a set of N variables, or ele- 
ments, which each have two possible values (say and 1). These elements 
interact with each other according to some given interaction rules, specified 
through a set of Boolean coupling functions that determine the variables at 
the next time-step, and thereby give the dynamics of the system. Such a dis- 
crete stepping of a set of Boolean variables, also known in general terms as 
a Boolean network, is of potential interest in several different fields, rang- 
ing from gene regulation and control, to modeling democracy and social 
organization, to understanding the behavior of glassy materials. 

The models were originally studied primarily for their biological interest, 
specifically by Stuart Kauffman who introduced the so-called N-K model 



in the context of gene expression and fitness landscapes in 1969 (Kauff- 
man [1969j , |1974| , |l995| , |l993j , |1990| , |l984| ). Since Kauffman's original work 



the scientific community has found a broad spectrum of applicability of 
these models. Specific biological problems studied include cell diffcrentia- 
tion ([Huang and Ingbcr [200(|), immune response ( Kauffman and Wcin 



berger [1989[), evolution flBoraholdt and Sneppcn [1998| ]; |Zawidzki [1998( | 
Bornholdt and Sncppcn [2000| ] ; Ito and Gunji [1994|), regulatory networks 
( Bornholdt and Rohlf [2000|] ) and neural netwo rks (Wang, Pichler, and Ross 



[1990|;|Dcrrida Gardner, and Zippclius [1987| ; [Kiirtcn [1988a{ ; |Bornholdt 



and Rohlf [200C| ] ) . In the first two examples, the basic binary element might 



be a chemical compound, while in the last it might be the state of firing of 
a neuron. A computer scientist might study a similar set of models, calling 
the basic elements gates, and be thinking about the logic of computer de- 
sign (Atlan, Fogclman-Soulic, Salomon, and Weisbuch [1981]; Lynch [1995[ ) 
or optimization ( Lee and Han [1998|; [Btauffcr [1994 ]). Earlier work in the 
mathematical literature ( Harris [1960| |; Metropolis and Ulam [1953 1) stud- 
ied random mapping models, which are a subset of the models introduced 
by Kauffman. This same kind of problem has also drawn consi derable at 



tention from physicists interested in the d evelopment of chaos (Glass and 



Hill [1998]; Luquc and Sole [1998, 1997a|; Kiirtcn and Beer ]1997|; Mestl, 


Baglcy, and Glass ]199 


7|; Baglcy and Glass [1996 1; I 


3hattacharjya and Liang 


[1996b]; Lynch [1995| 


) and also in problems associated with glassy and 


disordered materials ( 


Dcrrida and Flyvbjerg [198 


51; 


Dcrrida and Pomeau 


[1986|; Dcrrida [1987b 


; Dcrrida and Flyvbjerg [1987a 


]). In these examples, 



the Boolean element might be an atomic spin or the state of excitation of 
a molecule. 

In some sense, the type of Boolean networks introduced by Kauffman can 
be considered as a prototype of generic dynamical system, as they present 
chaotic as well as regular behavior and many other typical structures of 
dynamical systems. In the thermodynamic limit N — > oo, there can be 
"phase transitions" characterized by a critical line dividing chaotic from 
regular regions of state space. The study of the behavior of the system at 
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and near the phase transitions, which are attained by changing the model- 
parameters, has been a very major concern. 

As we shall describe in more detail below, these models are often studied 
in a version in which the couplings among the Boolean variables are picked 
randomly from some sort of ensemble. In fact, they are often called N-K 
models because each of the N elements composing the system, interact with 
exactly K others (randomly chosen). In addition, their coupling functions 
are usually picked at random from the space of all possible functions of K 
Boolean variables. Clearly this is a simplification of real systems as there 
is no particular problem which has such a generically chosen coupling. All 
real physical or biological problems have very specific couplings determined 
by the basic structure of the system in hand. However, in many cases the 
coupling structure of the system is very complex and completely unknown. 
In those cases the only option is to study the generic properties of generic 
couplings. One can then hope that the particular situation has as its most 
important properties ones which it shares with generic systems. 

Another simplification is the binary nature of the variables under study. 
Nevertheless, many systems have important changes in behavior when 
"threshold" values of the dynamical variables are reached (e.g. the synapses 
firing potential of a neuron, or the activation potential of a given chemical 
reaction in a metabolic network). In those cases, even though the variables 
may vary continuously, the binary approach is very suitable, representing 
the above-below threshold state of the variables. The Boolean case is par- 
ticularly favorable for the study of generic behavior. If one were to study 
a continuum, one would have to average the couplings over some rather 
complicated function space. For the Booleans, the function space is just a 
list of the different possible Boolean functions of Boolean variables. Since 
the space is enumerable, there is a very natural measure in the space. The 
averages needed for analytic work or simulations are direct and easy to 
define. 

In addition to its application, the study of generic systems is of mathe- 
matical interest in and for itself. 

1.1 Structure of Models 

Any model of a Boolean net starts from N elements {<ti, a%, ■ ■ ■ , &n}, 
each of which is a binary variable <7j G {0, 1}, i = 1, 2, . . . , N. In the time 
stepping, each of these Boolean elements is given by a function of the other 
elements. More precisely, the value of Oj at time t + 1 is determined by the 
value of its Ki controlling elements cr^u), < J j 2 (i)i ■•■ i <J j K (i) a ^ time t. In 
symbols, 

o-i (t+ 1) = fi(a- h (i) (i) , cTj- 2 (i) (t),..., a jK . (i) (i)) , (1.1) 

where fi is a Boolean function associated with the i th element that depends 
on Ki arguments. To establish completely the model it is necessary to 
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specify: 



the connectivity Ki of each element, namely, how many variables will 
influence the value of every o~i ; 

the linkages (or couplings) of each element, which is the particular set 
of variables Uj 1 u-s , <jj 2 (,-),..., Oj K u\ on which the element o~i depends, 
and 



the evolution rule of each element, which is the Boolean function fi 
determining the value of ai(t + 1) from the values of the linkages 
Cji^C*), °~j 2 (i){t), • ■ ■ , o- jK _ 

Once these quantities have been specified, equation ( p_.l| ) fully determines 
the dynamics of the system. In the most general case, the connectivities 
Ki may vary from one element to another. However, throughout this work 
we will consider only the case in which the connectivity is the same for all 
the nodes: Ki = K, i = 1, 2, . . . , TV. In doing so, it is possible to talk about 
the connectivity K of the whole system, which is an integer parameter by 
definition. It is worth mentioning though that when Ki varies from one 
element to another, the important parameter is the mean connectivity of 
the system, (K), defined as 

1 N 

i=i 

In this way, the mean connectivity might acquire non-integer values. Scale- 
free networks (Strogatz [2001]; Albert and Barabasi [200 if ), which have a 



very broad (power-law) distribution of Ki, can also be defined and charac- 
terized. 

Of fundamental importance is the way the linkages are assigned to the 
elements, as the dynamics of the system both qualitatively and quanti- 
tatively depend strongly on this assignment. Throughout this paper, we 
distinguish between two different kinds of assignment: In a lattice assign- 
ment all the bonds are arranged on some regular lattice. For example, the 
K control elements o-j 1 u\ , o-j 2 (j) , . . . , o~j K u-\ may be picked from among the 
2d nearest neighbors on a d dimensional hyper-cubic lattice. Alternatively, 
in a uniform assignment each and every element has an equal chance of 
appearing in this list. We shall call a Boolean system with such a uniform 
assignment a Kauffman net. (See Fig. |l.l| .) 

Of course, intermediate cases are possible, for example, one may consider 
systems with some linkages to far-away elements and others to neighboring 



elements. Small-world networks ( Strogatz [2001 1) are of this type 



For convenience, we will denote the whole set of Boolean elements {ai(t), 
02 (i), . . . , cr7v(i)} by the symbol E f : 

St = {<Ti(t) J o- 2 (t) I ...,ff JV (t)}; (1.2) 
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FIGURE 1.1. The different kinds of linkages in a one dimensional system, (a) In 
the Kauffman net the linkages of every element o"i are chosen at random among 
all the other elements ci ... on- (b) In a completely ordered lattice, the linkages 
are chosen according to the geometry of the space. In the case illustrated in this 
figure, Oi is linked to its first and second nearest neighbors. 



E t represents then the state of the system at time t. We can think of 
E t as an integer number which is the base-10 representation of the binary 
chain {<7i(i), 02 (i), . . . , <7jv(i)}. Since every variable ai has only two possible 
values, and 1, the number of all the possible configurations is ft = 2 , 
so that Et can be thought of as an integer satisfying < Et < 2 . This 
collection of integers is the base-10 representation of the state space of the 
system. Although it is not essential for the understanding of the underlying 
dynamics of the network, this integer representation proves to be very useful 
in the implementation of the computational algorithms used in numerical 
simulations (at least for small values of N). 

A note of caution is relevant at this point. We should distinguish the 
purely Boolean model described in this work from Kauffman's N-K land- 
scape model, which provides a description of fitness landscapes by including 
a fitness function to be optimized. We are not going to review on fitness 
landscapes since abundant literatur e already exists on this topic ( |Wilkc 
Ronncnwinkcl, and Martinctz [200i| ; |Kauffman [19951 , |l993| ). 
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TABLE 1.1. Illustration of two Boolean functions of three arguments. The first 
function is a particular random function, whereas the second one is a canalizing 
function of the first argument o\. When this argument is 0, the output of the 
function is always 1, while if <j\ = 1, the output can be either or 1. 



1.2 Coupling Functions 

The arguments of the coupling functions fi{o'j 1 ^, ■ ■ ■ , <7j K ^) can take on 
2 K different values. One specifies the functions by giving, for each of these 
values of the arguments, a value to the function. Therefore there are a total 
of 

H = 2 2 * (1.3) 



different possible functions. In Table 1.1 we give two examples of coupling 
functions for the case K = 3. There are 2 3 = 8 configurations of the ar- 
guments crj^i), (7 j 2 (i)i a j3(i)i an( l f° r cacn onc of these configurations, the 
function fi can acquire the values 1 or 0. For K = 3 there are 2 23 = 256 
tables similar to the one shown in Table 1, one for each Boolean coupling 
function. Different tables differ in their assignments of 0's and l's. If we as- 
sign a probability or weight to each of these functions, one gets an ensemble 
of possible couplings. 

Possible ensenble choices abound. One ensemble used extensively by 
Kauffman and others is the uniform distribution in which all functions 
are weighted equally. Alternatively, a magnetization bia^ may be applied 
by weighting the choice of functions with an outcome with a probabil- 



ity p, and the outcome 1 with a probability 1 — p (see, for example Bas- 
tolla and Parisi [1998b| ). Onc may also give different weights to particular 
types of functions. For example, one can consider only forcing functions 
or canalizing functions ( [Btauffcr [1987a ]; Kauffman [1969, p~984 |), in which 



the function's value is determined when just one of its arguments is given 



a specific value. The second function shown in Table 1.1 is a canalizing 



The word magnetization comes from the possibility of identifying each element with 
an atomic spin, which is a very small magnet. 
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TABLE 1.2. Boolean functions for K = 1. The first two functions form the class 
.4 of constant functions, Ao(u) = and Ai(cr) = 1. The other two functions form 
class B which consist in identity Bi(a) = a and negation Bn(c) = ^<r 
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TABLE 1.3. Boolean functions for the case K = 2. The 16 functions can be 
arranged in four different classes which differ in their symmetry properties (see 
text). 



function. Another possibility is to specify the value of the function in or- 



der to simulate the additive properties of neurons ( 


Bornholdt and Rohll| 


[2000 


|; iGenoud and Mctraux [199£]; Cheng and Tittcrington [1994 


1; Wang, 


Pichlcr, and Ross [199C|; 


Kiirten [1988a 


|; Dcrrida, Gardner, and Zippclius 


[1987|). 



Here we enumerate some of the coupling functions occurring for different 
values of the connectivity K. 

• For K = there are but two functions, corresponding to the two pos- 
sible values of a Boolean variables: tautology / = 1 and contradiction 
/ = 0. Together these two functions form a class which we might call 
A. 

• For K = 1, in addition to the class A, there exists another class B in 
which /(a) can take on the value a, called identity, and the value ->a, 
called negation. Thus there are a total of four functions, represented 



as columns in Table 1.2 



• The situation for K = 2 has been particularly carefully studied. Here 
there are four classes of functions /(cri,CT2) (Lynch [1993b|; Copper 



smith, Kadanoff, and Zhang [2001aQ. Each class is invariant under 
making the interchange «-► 1 in either arguments or value of /. The 
classes are A (two constant functions), B\ (four canalizing functions 
which depend on one argument), Bi (eight canalizing functions which 
depend on two arguments), C, (two no n-ca nalizing functions). These 



functions are explicitly shown in Table 1.3 



Several calculations have been done by giving different weights to the 
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different classes (see for example Lynch [1995[ ; Btauffer [1987a]) 



1.3 The Updates 

Once the linkages and the /j's are given, one is said to have defined a 
realization of the model. Given the realization, one can define a dynamics 
by using equation ( |1.]| ) to update all the elements at the same time. This 
is called a synchronous update. In this paper, we assume a synchronous 
update unless stated otherwise. Alternatively, one may have a serial model 
in which one updates only one element at a time. This element may be 
picked at random or by some predefined ordering scheme. 
Additional choices must be made. One can: 

1. Keep the same realization through all time. We then have a quenched 
model. 

2. Pick an entirely new realization after each time step. The model is 
then said to be annealed^ 

3. Employ a genetic algorithm in which the system slowly modifies its 



realization so as to approach a predefined goal (Bornholdt and Sncp- 
pen [20001 ; [Stern [1999j ; |Stauffer [1994| ). 



4. Intermediate choices are also possible ( [Baillic and Johnston [1994 ] 



Bastolla and Parisi [1996]) 



Almost always, we shall regard the real system as one which updates 
synchronously and which is quenched so that the interactions are fixed for 
all time. The annealed and sequential models will be regarded as approxi- 
mations which can provide clues to the behavior of this "real system" . The 
quenched model has time-independent dynamics describing motion within 
the state space of size f2 = 2 N . One iterates the model through time by 



using equation (1.1) and thereby obtains a dynamics for the system. Each 
of the Q different initial conditions will generate a motion, which will even- 
tually fall into a cyclical behavior. 



1.4 Symmetry Properties 

Typically, each individual realization of these models shows little or no 
symmetry. However, the average over realizations has quite a large symme- 
try group, and the symmetry is important to model behavior. For example, 



2 These terms have been borrowed from the physics of alloys in which something which 
is cooled quickly so that it cannot change its configuration is said be be quenched, and 
something which is held at a high temperature for a long time so that it can respond to 
its environment is described as annealed. Hence these terms are applied to situations in 
which one wants to distinguish between problems with fixed versus changing interactions. 
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the random mapping model ( Harris [1960( |), which is the K — > oo limit of 
the N-K model of the Kauffman net, has a full symmetry under the in- 
terchange of all states forming the state space. For finite values of K, the 
Kauffman net is symmetric under the interchange of any two basic elements. 
One can also have a symmetry under the interchange of the two values of 
each individual clement if one chooses the couplings at random, or with 
the ap propriate sym m etry. One can use dynamics that have rever sal sym- 
metry ( Harris [196C ]; |Coppcrsmith, Kadanoff, and Zhang [2001a , bf|), and 
that choice will have a profound effect upon the structure of the cycles. 



1.5 Outline of Paper 

To define fully the object of study, one must describe the dynamical process 
and the class of realizations that one wishes to study. For example, one can 
fix N and K, and study all of the properties of all realizations of frozen 
systems with those values. One might pick a realization at random among 
all possible linkages and functions, develop its properties. Then one would 
pick another realization, and study that. Many such steps would give us 
the properties of the frozen system averaged over realizations with a given 
N and K. What properties might we wish to study? 

In the next chapter, we describe the gross information transfer through 
the system by describing how the system will respond to a change in initial 
data or couplings. There are three different phases that have qualitatively 
different information transfer properties. We argue that the Kauffman net, 
which can transfer information from any element to any other, is qualita- 
tively different from lattice systems, in which the information transfer oc- 
curs through a d-dimensional space. We argue that the information transfer 
on the lattice is qualitatively, and even quantitatively, similar to the kind 
of information How studied in percolation problems. 

Chapter 3 is concerned with the temporal recurrence in the network as 
reflected in the statistical properties of its cycles. Here we base our argu- 
ments upon a discussion of two important limiting cases, K = 1, and very 
large K. The first case is dominated by situations in which there are a 
few short linkage loops. In the second, the Kauffman net shows a behavior 
which can be analyzed by comparison with a random walk through the 
state space. The distribution of cycle lengths is qualitatively different from 
any of the quantities that are commonly studied in percolation problems. 
So the results characterizing the cycles are substantially different from the 
behaviors usually studied in phase transition problems. We argue in ad- 
dition that the cycles of Kauffman nets and of networks on d-dimensional 
lattices differ substantially. 
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2 Information Flow 

2.1 Response to Changes 

The first thing to study in an N-K model is its response to changes. This 
response is important because the actual values of the elements often do 
not matter at all. If, for example, we pick the functions fi at random 
among the class of all Boolean functions of K variables, then the ensemble 
is invariant under flipping the value of the «th element. In that case, only 
changes matter, not values. 

In computer studies, such changes can be followed quite simply. One 
follows a few different time developments of systems that are identical 
except for a small number of selected changes in the coupling functions or 
initial data, and sees how the differences between the configurations change 
in time. One can do this for two configurations or for many, studying pair- 
wise differences between states, or things which remain identical across all 
the time-developments studied. 

2.1.1 Hamming Distance and Divergence of Orbits 

For simplicity, imagine starting out from two different possible initial states: 

E o = {0i(O),(7 3 (O),... ,^jv(0)} S = {&i(0),cr a (0),... ,a N (0)} (2.1) 

which differ in the values of a few elements. One can then observe the time- 
development of these configurations under the same dynamics, finding, for 
example, the distance D(t) between the configurations as a function of time 

N 

D(t) = Y^(*i(t)-ti(t)) 2 - (2-2) 

i =1 

If the transfer of information in the system is localized, this distance will 
never grow very large. If however, the system is sufficiently chaotic so that 
information may be transferred over the entire system, then in the limit of 
large N this Hamming distance can diverge for large times. 

Another interesting measure is the normalized overlap between configu- 
rations, a(t), defined as 

a(t) = 1 — N~ 1 D(t). (2.3) 

One will wish to know whether a goes to unity or a lesser value as t — > oo. If 
the overlap always goes to unity, independently of the starting states, then 
the system cannot retain a nonzero fraction of the information contained 
in its starting configuration. Alternatively, when a(oo) is less than unity, 
the system "remembers" a nonzero fraction of its input data. 
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2.1.2 Response to Damage 

So far, we have considered the system's response to changes of the initial 
data. One can also attack the quenched problem by considering two sys- 
tems, each with the same initial data, but with a few of the linkages varied. 
Then one can ask: given such "damage" to the system, how much do the 
subsequent states of the system vary? Do they become more and more alike 
or do they diverge? What is the likelihood of such a divergence? 

These considerations of robustness — both to damage and to variation in 
initial data — are very important for the evaluation of the effectiveness of a 
network, cither for computations or as part of a biological system. There 
have been fewer studies of the effect of damage than that of initial data. 



Usually the two types of robustness o ccur together ( Luque and Sole [200C ]; 



Dc Sales, Martins, and Stariolo [1997]). D amage has been studied for itself 



(Btauffcr [1994]; |Corstcn and Poole ftggg]) 



2.2 Percolation and Phase Behavior 
2.2.1 Percolation of Information 

In the limiting case in which ./V approaches infinity, the different types 
of N-K models all show three different kinds of phases, depending upon 
the form of information transfer in the system. If the time development 
transfers information to a number of elements that grows exponentially in 
time, the system is said to be in a chaotic phase. Typically, this behavior 
occurs for larger values of K, up to and including K = N. If, on the 
other hand, a change in the initial data typically propagates to only a 
finite number of other elements, the system is said to be in a frozen phase. 
This behavior will arise for smaller values of K, most especially K = 0, 
and usually K = 1. There is an intermediate situation in which information 
typically flows to more and more elements as time goes on, but this number 
increases only algebraically. This situation is described as a critical phase. 

When the linkages and the hopping among configurations are sufficiently 
random, one can easily perform a quite precise calculation of the boundary 
which separates these phases. Imagine starting with a state So, containing 
a very large number, TV, of Boolean elements, picked at random. Imagine 
further another configuration So in which the vast majority of the elements 
have the same value as in So, but nevertheless there are a large number of 
elements, picked at random, which arc different. The Hamming distance at 
time zero, D(Q), is the number of changed elements. Now take the simplest 
N-K system in which all the linkages and the couplings are picked at 
random. On average, a change in a single clement will change the argument 
of K functions, so there will be KD(0) functions affected. Each of these 
will have a probability one half of changing their value. (The functions after 
all are quite random.) Thus the Hamming distance after the first time step 
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will be 

D(l) = 0.5KD(0). 

If the couplings and connections are sufficiently random, then at the start of 
the next step, the newly changed elements and their couplings will remain 
quite random. Then the same sort of equation will apply in the next time 
step, and the next. Just so long as the fraction of changed elements remains 
small, and the randomness continues, the Hamming distance will continue 
to change by a factor of K/2 so that 

D{t + 1) = 0.5KD(t), 

which then has the solution 



D(t) = D(0)exp[iln(0.5iT)]. (2.4) 

For K > 2 the number of changed elements will grow exponentially, for 
K < 2 it will decay exponentially, and for K = 2 there will be neither 
exponential growth nor decay, and the behavior will be substantially influ- 
enced by fluctuations. Thus, by varying the value of the connectivity, the 
system sets down into one of the three following phases: 

• Chaotic (K > 2), the Hamming distance grows exponentially with 
time. 

• Frozen (K < 2), the Hamming distance decays exponentially with 
time. 

• Critical (K c = 2), the temporal evolution of the Hamming distance 
is determined mainly by fluctuations. 

In deriving equation (eT3) we have assumed that the coupling functions fi 
of the system acquire the values and 1 with the same probability p = 1/2. 
Nonetheless, as we will see below, the chaotic, frozen and critical phases 
are also present in the more general case in which the coupling functions fi 
evaluate to and 1 with probabilities p and 1 — p respectively. For a given 
value of p, there is a critical value K c (p) of the connectivity below which 
the system is in the frozen phase and above which the chaotic phase is 
attained. Conversely, for a given connectivity K > 2, a critical value p c (K) 
of the probability bias separates the chaotic and the frozen phases. 

The behavior of the normalized Hamming distance, D(t)/N, can be seen 



in figures 2.1 and 2.2, which respectively are for the Kauffman net and 
a two-dimensional lattice system. In both cases the system has N = 10 4 
elements and the connectivity is K = 4. The linkages of every element of 
the Kauffman net are chosen randomly, whereas in the two-dimensional 
lattice each element receives inputs from its four nearest neighbors. Both 
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FIGURE 2.1. Hamming distance as a function of time for a Kauffman net com- 
posed of N = 10000 elements and connectivity K = 4. (a) Log-log graph show- 
ing the Hamming distance for the three different regimes of the system: frozen 
(p = 0.05), critical (p c = (1 - \/T/2)/2 ~ 0.1464) and chaotic (p = 0.4). In all 
the cases the initial Hamming distance was D(0) — 100. (b) Hamming distance 
for the critical phase (p = p c ) but in a non- logarithmic graph. Note that the 
Hamming distance initially decreases, and then it rises again to saturate at a 
constant value that depends weakly on system size. 



figures contain three curves, with the parameter p picked to put the systems 
into the three diffe rent phases. For the Kauffman net p c = (1 — y/l/2)/2 
(see equation (2.8) below). The value of p c is not very well known for the 
two-dimensional lattice, but the best numerical estimations in d icate that 
p c ~ 0.29 for an infinite lattice (|Stauffer [1988|; |Stolzlc [1988| ; |Weisbuch 



and Stauffer [1987j ; |Stauffer [1987rj| ; [Dcrrida and Stauffcr [1986| ) . For finite 



lattices the value of p c has been defined as the average over realizations of 
the value of p at which a first cluster spanning the whole net appears (Lam 
[1988| ). For a 100 x 100 lattice this value is p c ~ 0.27. 



In the frozen phase the distance has an initial transient but then quite 
quickly approaches an asymptotic value. In the chaotic phase the distance 
shows an exponential rise followed by a similar saturation. These behaviors 
are almost the same for both the Kauffman net and the two-dimensional 
lattice. On the other hand, in the critical phase the behavior of the Ham- 
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FIGURE 2.2. Hamming distance in a two-dimensional lattice composed of 
N = 100 x 100 elements. Every node in the lattice receive inputs from its four 
first nearest neighbors (K = 4). The three curves (in log-log scale) show the 
behavior of the Hamming distance in the three regimes: frozen (p = 0.1), critical 
(p c = 0.27) and chaotic (p — 0.4). Note that in the critical phase the Hamming 
distance initially increases algebraically and then saturates at a constant value 
that depends on system size. 



ming distance is very different in these two systems. In the Kauffman net 
the distance initially decreases and then increases again, asymptotically ap- 
proaching a constant value that depends weakly on system size. In contrast, 
for the lattice the Hamming distance initially grows and then saturates. We 
will see later that, within the framework of the annealed approximation, 
the normalized Hamming distance for an infinite Kauffman net approaches 
zero monotonically in both the frozen and the critical phases (exponen- 
tially in the frozen phase, and algebraically in the critical phase). As far as 
we can tell, the non-monotonic behavior of the Hamming distance in finite 



systems at K c shown in Fig. 2.1b has not yet been explained 



2.2.2 Limitations on the mean field calculation 



Let us look back at the argument which led to equation (2.4). Calculations 
like this, in which actual system properties are replaced by average system 
properties are in general called "mean field" calculations. 



Naturally, the results derived as equation (2.4) depend crucially upon the 
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assumptions made. The derivation follows from the assumption that the fa' a 
in each step are effectively random. (See also the derivations of equations 
( f2~6| ) and ( 2.7) below, which also depend upon the randomness assumption.) 
The randomness will certainly be correct in the annealed situation, in which 
the couplings are reshuffled in each step. It will also be true in Kauffman 
net in the limiting situation in which K = oo. In that case, information 
is spread out over the entire system and thus has a very small chance of 
correlation among the /j's. The Kauffman net has a likely configuration that 
permits the replacement of the actual values of the fa's by their statistical 
distribution ( Hilhorst and Nijmcijer [1987 1 ) . However, the approximations 
used here will not always work. Specifically, they fail in all kinds of finite 
N situations, or in situations in which the linkages are arranged in a finite- 
dimensional lattice. In that case, the assumed randomness of the fa does not 
hold, because their arguments are not random, and the derived equations 
will not work. To see this in the simplest example choose K = N = 1 with 
quenched couplings. A brief calculation shows that for any one of the four 
possible fa's, after a one-step initial transient a(t + 2) = a(t). That does not 
agree with equation ( |2.6| ) derived below. In fact, for any finite dimension 
and linkages which involve short-range couplings, the overlap is not unity 
at long times even in the frozen phase. More generally, if the system is 
put onto a finite dimensional lattice, or if the functions are not picked 
at random, or if the initial elements are not random, couplings initially 
used can be correlated with couplings used later on. Then the information 
transfer will be different and equation (2.4) will fail. 

However, the principle that there can be a phase transition in the kind 
of information flow remains quite true for a d-dimensional lattice, and for 
other ensembles of coupling functions. 



2.2.3 Connections to percolation problems 

The transfer of information just described is quite similar to the trans- 
fer which occurs in a typical phase transition problem. Gcneri cally, thes e 



problems have three phases: ordered, critical, and disordered ( Ma [1976 



Kadanoff [2000 1 ) . The bare bones of s uch an informa tion transfer problem 



is described as a percolation problem ( Btauffer [1985 1 ) 



In one sort of percolation problem one considers a lattice. Each bond 
or link of the lattice is picked to be either connected or disconnected. The 
choice is made at random and the probability of connecting a given link 
is picked to be q. Now one asks about the size of the structure obtained 
by considering sets of links all connected with one another. For small val- 
ues of q, these sets tend to be small and isolated. As q is increased, the 
clusters tend to get larger. At sufficiently large values of q, one or more 
connected clusters may span the entire lattice. There is a critical value of 
the probability, denoted as q c , at which such a spanning cluster just begins 
to form. Properties of the large clusters formed near that critical point 
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have been studied extensively ( Stauffcr [1985[ ). The resulting behavior is 
"universal" in that for an isotropic situation, the critical properties depend 
only upon the dimensionality d of the lattice, at least when d is less than 
four. For d > 4, the percolating system obeys mean field equations. When 
the symmetry properties change, the critical behavior does change. For 
example, a system with directed percolation (Stauffcr [1985 1 ; pwezarck^ 



Rechnitzer, and Guttmann [1997 1) has the information flow move preferen- 



tially in a particular direction. The critical behavior of directed percolation 
is different from that of ordinary percolation. It is attractive to explore the 
connection between the phase transition for percolation, and the one for 
the N-K model. 

Several authors have constructed this connection in detail. For example, 



Hansen [1988a | looked at the N-K model on lattices for 2, 3 and 4 dimen- 



sions and demonstrated numerically that the phase transition occurred 
when the information transition probability reached the critical value for 
percolation on the corresponding lattice. |5tolzle [1988 1 showed the con- 
nection to percolation for both sequential and parallel updating for a two 
dimensional lattice. However, Hansen [1988b[ took very special forms of 
the connections, using only rules with connectivities of the form of a logi- 
cal "or" . This did not result in behavior like simple percolation but instead 
quite a different phase transition problem, related to the behavior of diodes 
and resistors. At roughly the same time, Stauffcr [1988] ] indicated a close 
numerical correspondence to the two dimensional percolation problem both 
in p c and also in the fractal dimension of the spanning cluster. (For p c see 
also Lam [1986| ].) He also defined an order parameter, essentially a Ham- 
ming distance, that, when plotted as a function of (p — p c ), looked just 
like a typical critical phenomena result. He argued that the N-K model is 



equivalent to directed percolation. More specifically, Obukhov and Stauf- 



fcr [1989 1 argued that the quenc hed problem has a critical point which is 
in the same "universality class" ( Kadanoff [2000| ) as directed percolation. 
This would imply that the critical point is essentially the same as that of 
the directed percolation problem. The qualitative properties of both the 
ordered and the frozen phases would also be essentially similar in the per- 
colation case and the N-K lattice. In the latter situation, the preferred 
direction would be that of the "time" axis. The structure of both would 
vary with dimensionality and become like that of mean field theory above 
four dimensions. This is the same mean field theory which describes the 
Kauffman net. Thus, the behavior of information transfer in N-K problems 
was mostly understood in terms of percolation. 



2.3 Lattice versus Kauffman net 

We can now point to an important difference between systems in which 
all elements are coupled to all others, as in the Kauffman net, and lattice 
systems in which the elements which are "close" to one another are likely 



Version Jan. 21, 2002 



Edited by : Typeset on 5 February 2008 - 13h50 



1. Dynamics in Boolean Networks 17 



to be mutually coupled. "Closeness" is a reciprocal relation. If a is close 
to b, then b is also close to a. Therefore, close elements are likely to be 
coupled to one another and thereby form a closed linkage loop. Any large- 
TV lattice system might be expected to form many such loops. When K is 
small, different spatial regions tend to be unconnected and so many different 
modules will formP]. The dynamics of the elements in different modules are 
independent. In contrast, in a Kauffman net, influence is not a reciprocal 
relation. If element o~j appears in the coupling function fi associated with 
element Cj, there is only a small chance, proportional to (K/N), that <7j will 
appear in fj. For large N and small K, the probability that a given element 
will participate in a linkage loop will be quite small, so there will then be 
a small number of modules. When K is small, the number of modules in 
uniformly coupled systems grows more slowly than the system size, while 
in lattice systems the number of modules is proportional to the size of the 
system. This distinction will make a large difference in the cycle structure.^ 
For the flow of information the difference between the two kinds of nets is 
more quantitative than qualitative. One can see the similarity between 
them by comparing the curves shown in figures 2.1 and 2.2. 



2.4 Calculations of overlap and divergence 

Before coming to a careful description of the phases, we should describe 
more fully the kinds of calculation of overlap that can be performed. Equa- 
tion ( p,. 4) is just the beginning of what can be done with the trajectories 
of states in this problem. In fact, exactly the same logic which leads to 
that equation can give a much more general result. If the overlap between 
two states at time t is a(t), and if the elements which are different arise at 
random, then the probability that the arguments of the function fi will be 
the same for the two configurations is 

p=[a(t)] K . (2.5) 

If all arguments are the same, then the contribution to the overlap at time 
t + 1 is 1/N. (The N arises from the normalization of the overlap.) If 
one or more arguments of the coupling function are different in the two 
configurations, and the functions fi are picked at random, then the chance 



3 A module is a loop or more complex topology of dependencies in which all the 
functions a rc n on-constant, plus all other elements that are influenced by that structure. 



See section 3.1 . 

4 The interested reader will recall that in quantum field theory and statistical me- 
chanics, mean field theory appears in a situation in which fluctuations are relatively 
unimportant. This will arise when the analog of linkage loops make a small contribu- 
tion. Then, the physicist puts the effect of loops back into the problem by doing what 
is called a loop expansion. Since they both expand in loops, the percolation mean field 
theory and the mean field theory of the Kauffman net are essentially the same. 
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of having the same functional output is 1/2 and the contribution to the 
overlap is 1/(27V). Since there are N of such contributions, weighted by 
p and 1 — p respectively, the equation for the overlap defined by equation 
© is 

a(t + l)= [l + [a(t)] K ]/2. (2.6) 

There are several possible variants of this equation. For example, if the 
different outcomes and 1 of the function fa are weighted with probabilities 
p and 1 — p respectively, to produce a sort of magnetization bias, then 
equation (f2~(i|) is replaced by ( Dcrrida and Pomcau [1986 1; Derrida and 



Stauffer [1986}] ) 

a(t+l) = l- [l-[a(t)] K ]/K c , (2.7) 

where K c is given in terms of p as 

K D =l/[2p(l-p)]. (2.8) 

In the limit t — > oo, a(t) asymptotically approaches the fixed point a*, 
which obeys, from equation (|2.7D 



1- [l-[a*] K ]/K c . (2.9) 



We might expect equation (2.9) to reflect the three-phase structure of the 



problem, and indeed it does. Fig. 2.3a shows the graph of the mapping 
F(a) = 1 — [l — a K ] /K c for different values of K, and Fig. |2.3[ j shows the 



bifurcation diagram of equation (2.9). Both graphs were calculated with p 



chosen so that K c = 3. As can be seen, if K < K c there is only one fixed 
point a* = 1, whereas for K > K c the fixed point a* = 1 becomes unstable 
as another stable fixed point, a* ^ 1, appears. The value of the infinite 
time overlap a* describes the fraction of elements whose value is sensitive 
to the cycle entered. 

When K > K c . the system is chaotic; a* is less than one even when the 
starting points are quite close to one another. This reflects the fact that the 
system has a long-term behavior which can include some cyclic oscillations, 
the initial data determining the cycle entered. (We discuss the cycles in 
detail in the next chapter.) As K approaches K c from above, a* increases 
since fewer elements have final values which are sensitive to the initial data. 
On the other hand, for K < K c , the infinite time overlap is exactly one 
and therefore the proportion of elements whose final value depends upon 
the starting point is precisely zero. Thus, independently of the starting 
point, the system is always stuck in essentially the same configuration. 
This surprising result pertains to the Kauffman net. In contrast, for any 
lattice system with K > the final overlap is less than unity, reflecting the 
fact that the system can store information about initial data in some finite 
proportion of its elements. Such storage is impossible for the Kauffman net. 

Version Jan. 21, 2002 Edited by : Typeset on 5 February 2008 - 13h50 



1. Dynamics in Boolean Networks 19 



a(t+l) 



1 

0.9 
0.8 
0.7 
0.6 
0.5 



K=l 


K = 3 




' A" = 7 


i 






\ * 



0.2 0.4 0.6 

(a) 



0.8 



aft) 



1.00 



0.80 



0.60 




10 

(b) 



12 14 



16 



18 20 



FIGURE 2.3. (a) The mapping F(a) = 1 - [l - a K \ /K c (see Eq. (JfLty) for 
iv c = 3 and three different values of K (solid curves), corresponding to the 
three different phases of the syst em. The dotted line is the identity mapping, (b) 
Bifurcation diagram of equation ( |2.9| ). For K < K c the only fixed point is a* = 1. 
For K > K c the previous fixed point becomes unstable and another stable fixed 
point a* a appears. 



It is worth noticing that for a given value of K, the system can be put in 
one of the three different phases by adjusting the value of the proba bilit y 
bias p. The critical value of p is then obtained by inverting equation (2. 



which gives the critical line separating the frozen phase from the chaotic 



phase in the p~K parameter space, as illustrated in Fig. 2.4 



The rigidity of the Kauffman net was emphasized by Flyvbjerg (Fly- 



vbjerg [1988 , 198S |), who wished to establish the equation for the stable 



core, the set of variables whose values at time infinity do not depend upon 
initial data (see also Lynch [1995 , 1993a, EJ; Bastolla and Parisi [19981 1.) 
He calculated the time dependence of the proportion of variables which 
had already settled down, and thereby found a closed form for the size 
of the core. He found that in the ordered state of the Kauffman net the 
fraction of elements in the core is unity. Another kind of rigidity is studied 
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FIGURE 2.4. Phase diagram for the N-K model. The shaded area corresponds 
to the frozen phase, whereas the upper region corresponds to the chaotic phase. 
The curve separating both regions is the critical phase K c = [2p(l — p)]^ 1 ■ 



by counting the weak elements. In a particular realization, elements arc 
termed weak if changing their value never affects the long-term behavior 
of the system. Man y authors have studied these elements (see for example 
Lynch [1995 , 1993a ,B]). In the Kauffman net, but not in the lattice system, 



the proportion of weak elements is unity throughout the frozen phase. 
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3 Cycle Behavior 

For quenched systems the evolution functions fi are independent of time. 
For any initial condition, the system eventually returns to a previously 
visited point and then cycles repeatedly. The state space, which consists of 
the fl = 2 N configurations of the system, breaks up into different cycles and 
their basins of attraction, as shown schematically in Fig. |3.1[ Each initial 
condition will eventually fall in one of these cycles, which are determined 
by the evolution functions and the linkages. 




FIGURE 3.1. Schematic representation of the state space of the N-K model. 
Each state is represented as a bold point. Under the quenched dynamics, the state 
space is broken down into several cycles, or attractors, represented as circles. Each 
initial state eventually will end up in one of these cycles (the arrows show the 
direction of the flow) . The totality of states which evolve towards a given cycle, is 
the basin of attraction of that cycle. Note that there can be attractors consisting 
of only one point, corresponding to a cycle of period 1. 

The description of a cycle is, in some sense, much more complex than 
the description of orbit separation. In separation, one is dealing with a very 
gross property of the system: Different parts of it can behave independently 
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and additively in the separation process. We utilized this fact in our calcu- 
lations of overlaps and Hamming distances carried out in sections 2.2.1 and 
2.4 above. On the other hand, to get a cycle to close, each and every element 
must simultaneously return to its value at a previous time-step. This is a 
much more delicate process and makes for much more delicate calculations. 
As we saw, information flow in N-K systems is closely analogous with the 
well-studied behavior of percolation problems. In contrast, the behavior of 
cycles in N-K models is something special, not closely analogous to any of 
the usual calculations for the usual systems of statistical mechanics. When 
the N-K system forms cycles, one can ask many quantitative questions. 
One might wish to know the number of steps the system takes before it 
falls into a cycle (called the transient time) , and about the length of the 
cycle. For each cycle, there is a basin of attraction, which is the set of ini- 
tial states which eventually fall into that cycle. Therefore, one can develop 
a quantitative description of the distribution of cycle lengths, transient 
lengths, basin sizes, etc. either within a given realization or averaged over 
realizations. 



3.1 Linkage Loops and Cycles 

The structure of the linkages is an important ingredient in determining the 
nature of the cycles in the system. Each element can appear in the coupling 
functions for other elements. These in turn can appear in the couplings of 
others. For each element one can trace out its descendants, i.e., the elements 
it affects. Similarly, one can chain backward and find the ancestors for each 
element, i.e., the elements which affect it. For a cycle of length longer than 
one to exist, at least one element must be its own ancestor, and thus its 
own descendant. (If no such clement existed, one could trace back and find 
elements with no ancestors. They would then have fixed values and be sta- 
ble elements. The elements controlled by only them would also be stable. 
The line of argument would go forward until it was found that all elements 



were stable.) Fig. 3/2 illustrates the idea of ancestors and descendants for 
the case K = 1. As can be seen, the connections between different elements 
of the network gi ve rise to linkage loops and linkage trees, each tree being 



rooted in a loop. Flyvbjerg and Kjacr [1988 ] stress the importance of el 



ements which are their own ancestors by pointing out that any unstable 
element must be influenced by at least one such element. Of particular rel- 
evance are those elements which belong to a loop in which there are only 
non-constant functions. These self-influencing elements, and the unstable 
elements that they influence, are called the relevant elements of the sys- 



tem. They naturally form themselves into groups called modules (Thicffry 
and Romero [1999| ; |Bastolla and Parisi [1998a] ; ^awidzki [199^1 ) . Different 



modules do not influence one another and fully determine the cycle struc- 
ture. Notice that the chain of linkages by which an element can influence 
itself, namely its linkage loop, is completely defined by the linkages, i.e. 
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FIGURE 3.2. Linkage loops for an N-K net of N — 20 elements and connec- 
tivity K — 1. (a) Particular realization of linkage assignment in the net. The 
first column is the list of the N elements of the net, {cti, ai, . . . , crjv}, and the 
second column shows the particular linkage every element has been assigned, 
(b) Schematic representation of the linkage loops. Each arrow points towards 
the descendants of a given element. In this particular realization there are three 
modules, one of which consits of only one element, 02- 

by the specific assignment of the different elements in the fa's. Linkage 
loops do not describe the functions themselves. Consequently, such loops 
are necessary but not sufficient for the existence of non-trivial cycles. Only 
loops formed by relevant elements (with non-constant functions) are impor- 
tant in determining the properties of the limit cycles of the network. The 
number of such cycles and their leng ths depend crucially on the modular 
organization of the relevant variables ( Flyvbjerg and Kjacr [1988 ] ; Bastolla 
and Parisi [1998b , af). Consequently, linkage loops and relevant elements 



will form an important part of our further discussions. 



3.2 Phase Transitions for Cycles. 

In the previous chapter, we saw that N-K systems fall into different phases, 
depending upon how effectively they transfer information. The cycles can 
be quite different in the different phases. The chaotic phase is characterized 
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by very long cycles, in which the typical cycle- lengths grow as a power of 
the size, f2 = 2^, of the entire state space. Each cycle includes the par- 
ticipation of many different elements which undergo a strongly correlated 
dynamics. In this phase, the transients are similarly long and complex. In 
contrast, the frozen phase tends to have much shorter and simpler cycles, 
and also shorter transients in which the individual modules do not grow 
as the system gets larger. If the system sits on a lattice of low dimension, 
the different modules are spatially localized. Each module shows a strong 
correlation in its dynamics, but different modules move independently. The 
critical phase has larger modules and has a behavior intermediate between 
the frozen and the chaotic phase. These behaviors are known in consider- 
able detail. 

The remainder of this chapter describes in detail the behaviors of the 
different phases. We start by discussing the exactly solvable models which 
give a solid description of the limiting behavior of the phases. We then 
describe how information obtained from simulations and the exact solutions 
can be pieced together to give a qualitative description of the phases. 

3.3 Soluble Models. 

3.3.1 Independent Subsystems 

Here we do a preliminary calculation which will be of use in interpreting 
results involving cycles. Imagine a system composed of N independent sub- 
systems. Each subsystem has a probability pi of having a cycle of length /. 
We imagine that pi gets quite small for large I and ask what is the chance 
of finding a long cycle in the entire system. Notice that the chance of not 
finding a piece with a cycle of length I in the entire system is 

qi = (1 - Pl ) N « exp(-Npi). 

If then, pi varies exponentially with I (the justification for this assumption 
will be given below), namely, if 

pi = Aexp(-al), (3-1) 

then we might expect to find parts with all cycle lengths up to 

l mx = (In AO/a (3.2) 

(so that qi is not that small, say of order qi ~ e~ A ). To make a long cycle in 
the entire system, one puts together many sub-cycles of different lengths, 
li, i = 1, 2, . . . , N. The total cycle length, L, is the smallest number divisible 
by each of the ij's. Then L will be a product of all prime numbers, p r , which 
are less than l mxi each raised to a power s r which is the largest integer for 
which the inequality 

s7~ — 
.Pr . 
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is satisfied ([x] being the integer part of x). Hence, to a decent approxima- 
tion, the largest cycle length L mx will be 

L ps (I 

where tt(1) is the number of primes less than I, which can be estimated in 
the asymptotic limit of large I as |^| 

?r(0 = l/lnl. 

In the end then, the longest cycle length L mx obeys 
so that 

L mx « A^ 1/q . (3.3) 

We have reach the remarkable conclusion that even though the probabil- 
ity of long cycles in each component of the system falls exponentially, the 
typical maximum cycle length in the entire system depends algebraically 
upon the size of the system. This calculation does not apply directly to 
Kauffman nets because we have not accounted for the fact that the differ- 
ent modules have different distributions pi, but nonetheless it is instructive. 



3.3.2 K = 

The case in which K is zero is simple and uninteresting. After the first step, 
each element has a value which is completely determined by its function . 
Each element remains with the value it had at time one for all subsequent 
times. Thus the system is completely frozen. 



3.3.3 K = l 



In an important paper, Flyvbjerg and Kjacr [198§| ] analyze the structure of 



the case in which each coupling function depends upon the value of just one 
element. To analyze this case the authors focus upon those spins which are 
ancestors of themselves. As we have seen (see Fig. |3.2| ), each such element 
forms part of a coupling loop of a length which we will denote as to. No 
information may pass into such a loop from other parts of the system. Each 
element in the loop may affect others, but the affected elements either 
are constant or they inherit the cycle length of the coupling loop. In a 
lattice system, each element is coupled with a neighboring element. In any 
finite number of dimensions, there is a nonzero probability that two nearest 



5 A better approximation is tt(1) = l/(lnZ — 1) (Rosser and Schoenfeld [1962 1 ) . Also, 
in this reference it is shown that l/lnl < w(l) < 1.0423i/lnZ for all I > 10. 
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neighbor sites are inputs to each other. Therefore, a system with N sites 
will contain a number of loops that is proportional to N. In contrast, on the 
Kauffman net, the couplings are not to neighbors but randomly chosen from 
the whole system. When K = 1, the probability that two sites are inputs to 
each other is proportional to 1/-/V, and the average number of loops grows 



logarithmically with N (Bastolla and Parisi [1998a]). Flyvbjerg and Kjaer 



calculate the probability of observing n m loops of length m in a system of 
N — ► oo elements. Let rriT be the total number of elements contained in 
loops 

oo 

m,T = ^2 mn m . 

m— 1 

Then the distribution of {m, 112, ■ ■ ■} takes the form 

P( ni ,n 2 ,...) = ^eM-m 2 T /(2N)] f[ OJl^l, (3.4) 

m— 1 

Thus one can have reasonably long loops with a maximum likely length 
of order N 1 / 2 . The linkage loops do not determine the cycle structure. 
To know the number and kind of cycles, one has to know the coupling 
functions. For a given loop to make a non-trivial cycle, all the functions on 
it must be either identity or negation, for if there is a constant function in 
the loop, one element has a fixed value and it passes on its constancy to 
all the other elements in the loop. So think of a specific case. Let us have 
a loop, with m — 3, in which element 1 is the ancestor of clement 2 which 
is the ancestor of element 3 which is the ancestor of element 1. Let all the 
coupling functions be identity. Then the initial data just cycles around the 
loop. If the initial data is (ABC) then next is (CBA), etc. There are two 
cycles of length one in which all elements are identical and two of length 
three in which they are not. On the other hand, if all the coupling functions 
of this three elements are negation, there is one cycle of length two and one 
of length six. All the different situations are similarly easy to analyze. 

Despite the fact that the Kauffman net can contain very large linkage 



loops for K = 1, the cycles are reasonably short. Fig. 3.3a shows the distri 



bution o f cy cle lengths in a Kauffman net composed of TV = 15 elements, 



and Fig. 3.3 b shows the probability W(n) that an arbitrary state belongs to 
a basin of attraction of size n. For N — 15, there are 2 15 = 32768 different 
states. Nonetheless, as can be seen from the figure, the probability of find- 
ing a cycle with length bigger than 10 is negligibly small. The reason is that 
the probability that a loop is relevant (has only non-constant functions) de- 
creases exponentially as the loop length increases. In the usual calculation 
one assigns the K = 1 functions with equal weight. If either of the two 
constant functions are present in the loop, the cycle length is unity. If a is 
the probability of assignment of the constant functions and b = 1 — a is the 
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FIGURE 3.3. (a) Distribution of the probability that a cycle has length I, P(l), 
as a function of I for a Kauffman net of N = 15 elements and K = 1. Note that 
the probability of having cycles with lengths larger than 10 is rather small, (b) 
Probability P(w) that an arbitrary state belongs to a basin of attraction of size 
w, also for N = 15 and K = 1. 




probability of assignment of the other two functions, identity and negation, 
the probability of finding a cycle of length I is proportional to b l . In most 
calculations b is 1/2, so that the probability of really long cycles falls off 



exponentially in I (this is the justification of equation (3.1)). Flyvbjerg and 



Kjaer [198S | point out that in the special case with 6=1, the probability 



of long cycles falls algebraically. They speculate that the behavior in this 
limit might be, first, analyzable, and second, very similar to the behavior 
of the Kauffman net in the critical case (K = K c ) described below. As far 
as we know, these speculations remain unproven. 



3.3.4 K > N 

Another case in which it is possible to analyze the structure of the cycles 
in great detail is the one in which each coupling function depends upon all 
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the values of all the different elements in the system. [] This case has the 
simplifying characteristic that a change in a single spin changes the input 
of every coupling function. Therefore, one can analyze some features of the 
behavior of the system as if the system were annealed rather than quenched. 
In particular, one can calculate probabilities for hopping from configuration 
to configuration as if the system were undergoing a random walk through a 
space of size fi. The quenched nature of the system only asserts itself when 
the hopping takes onto a configuration previously visited. After that, one 
can be sure that the subsequent behavior will be cyclic. There are almost 
classical mathema tical analyses of this situation ( |Harris [196C ] ; Dcrrida and 



Flyvbjcrg [1987b | ) . We describe this case by considering the calculation of 
typical distributions of cycle lengths and of transients. Imagine a Kauffman 
net with K > N . Imagine that we start from a random configuration Eo 
at time t = 0. At subsequent times, we step by step follow the dynamics 
and essentially go through a random walk So, Si, £2,..., through the 
configuration space, which has size Q. This walk continues until we land 
upon a point previously visited. Let us define qt as the probability that 
the trajectory remains unclosed after t steps. If the trajectory is still open 
at time t, we have already visited t + 1 different sites (including the sites 
£0 and £*)■ Therefore, there are t+1 ways of terminating the walk at the 
next time step and a relative probability of termination p t = (t + 1)/Q. 
The probability of still having an open trajectory after t + 1 steps is 

Qt+l = - Pt) = <h ( 1 - -ff^j with 1o = 1, 

while the probability pt+i of terminating the excursion at time t + 1 is 

t + 1 

Pt+1 = ~cT qt ' 

To obtain P(L), the probability that a given starting point is in the basin 
of attraction of a cycle of length L, we note that a closure event at time t 
yields with equal probability all cycle lengths up to t. Therefore, 



t=L 

which, in the limit of large Q can be approximated by 

poo -, 

P(L)„ / JLe-xfc-iVpn)^ (3. 6 ) 



6 The resu l ts of this section have been kn own for quite some time; see references 
Harris [196C|; Dcrrida and Flyvbjcrg [1987b[ , which study what is called the random 



map model. Here one studies systems in which one has a random map from point to 
point in configuration space. The form of argumentation in this section closely follows 
Coppersmith, Kadanoff, and Zhang 
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It is also important to consider the probability P(m, L) of finding a cycle 
of length L after having gone through a precursor of length m, given by 



P(m,L) = ^q. 



m+L-l- 



(The factor 1/Q comes from the fact that only one point, £ m , of the state 
space can split the entire sequence {So, Si, . . . , S m , £ m +i, . . . , of 
m + L states into two pieces of lengths m and L respectively) In the limit 
of large f2, the previous expression can be approximated by 



P{m,L) 



exp[-(m + L) 2 /(2fi)] 



(3.7) 



The most important characteristic of the results ([3J]) and (^J) is that 
the typical cycle length and the typical precursor length are each of order 
VL 1 / 2 . Thus, a very small fraction of the total configurations participate 
in each transient or cycle, but nonetheless the cycles and the transients 
may be very long. There is another, and very nice, interpretation of the 
results just calculated. If the precursor length is zero, equation (3.7) gives 
the probability that our system will contain a cycle element in a cycle of 
length L. Since there are possible starting points, the average number of 
cycles of length L in our system is 



(N C (L)) 



exp[-L 2 /(2fi)] 



(3. 



Here the (• • • ) represents an average over realizations. An integration then 
gives us the information that the average total number of cycles is propor- 
tional to InO, or more precisely 7 



(N C ) = ^N + 0(1). 



(3.9) 



Fig. ^a shows the distribution P(L) for a Kauffman net with N = 
K = 10, while Fig. 3.4b shows the probability P(w) that a given starting 
point belongs to a a basin of attraction of size w. Note that most of the 
distribution P(w) is concentrated around w = 1 and w — Q, being nearly 
constant (zero) in between. The above reflects the fact that there are large 
fluctuations in the number of cycles in different realizations. If the num- 
ber of cycles were a typical extensive quantity, its median would be the 
same as its average value. Instead, for large N its median value is half the 
average, indicating that, in the Kauffman net, the average is dominated 
by a few situations with anomalously many cycles ( Harris [1960 1; Derrida 



3 Another result, often reported i n the literature (Zoli, Guidolin, Fuxe, and Agnat 
, or table 1 in Kauffman [1990 1 ) is (N c ) = N/e. We do not know the justification 



[l99flt 



for this, and suspect that it is wrong. 
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FIGURE 3.4. (a) Probability distribution of cycle lengths for a Kauffman net 
with N = K = 10. The solid curve is the theoretical result (3.5) and the broken 
line is the result of our numerical simulation, (b) Probability that a given starting 
point Eo belongs to a basin of attraction of size w, also for the case N = K = 10. 

and Flyvbjerg [1987b]). The typical size of the largest basin is given by 
dividing the entire space (of size O = 2 N ), by the number of cycles (of 
order N, which is actually of order one, given the variance involved with 
the various fluctuations). To form a basin of size ft = 2 N , one starts from 
the cycle with its typical size 2 JV / 2 and count backward finding the set of all 
first order ancestors, then backwards again to find their ancestors, which 
are second order ancestors of the cycle, and so on. Since these classes are 
non-overlapping and bounded in size, after some number of steps (perhaps 
of order 2 Af//2 , the typical precursor size) the classes will shrink in size and 
eventually one will terminate on one or a few of the most remote ancestors. 
The total basin size is the sum of the number of elements in these ances- 
tor classes. Since each configuration may have no ancestors, one ancestor 
or many ancestors, this backward-stepping, ancestor-counting process will 
have a character similar to a multiplicative random walk. As such it is an 
inherently highly fluctuating process. So we might expect huge fluctuations 



in the basin sizes (see Fig. 3.4 b) 
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3.4 Different Phases - Different Cycles 
3.4.1 Frozen Phases 

In the frozen phase, information typically propagates from a given element 
to only a very few other elements. Thus a change in initial data will typically 
affect the subsequent behavior of only a few elements. Similarly, damage to 
a single coupling function will produce changes which propagate to only a 
limited number of elements in the system. In these respects the lattice N-K 
models and the Kauffman net are very similar. However, in other respects 
they are very different. For example, in the lattice system, the number of 
relevant modules is proportional to N, so each initially different element 
has a nonzero chance of being in a relevant module, and each difference 
at time zero has a nonzero probability of developing into a difference at 
infinite time. In contrast, in the Kauffman net there are typically only 
a few short cycles and the overlap rapidly approaches unity, no matter 
what its initial value might have been. For the N-K model on the lattice 
there are many loops. In one dimension with nearest neighbor couplings 
all loops have length unity, but with longer range couplings and in more 
dimensions some loops might be quite long. Despite the fact that at any 
given point in the lattice it is exponentially improbable to find long loops, 



the argument of section 3.3.1 indicates that the average loop size might 



well vary algebraically with the size of the system. 

1. Lattice system. For small K, the lattice system falls into a phase 
in which there are many short cycles. The number of cycles grows 
as a power fi of the volume of the state space, f2; thus, the typical 
basin size will be The typical cycle length grows as a power 
of N, with multiplicative fluctuations of order IniV. The growth in 
cycle length occurs because different modules will have sub-cycles of 
different lengths. The entire cycle length is the product of the prime 
factors coming from the modules. Large primes arc exponentially un- 
likely in a given region, but the number of regions observed grows 
with N. 

Compare the time development coming from two starting configura- 
tions in which some small fraction of elements, spread out through 
the system, are different. These differences will each have a finite 
probability of causing a different cyclical behavior, localized in its 
particular region of the lattice. Independent of the size of the sys- 
tem, at long times the Hamming distance will go to a constant D*, 
proportional to -D(O) with a constant of proportionality of order one. 

2. Kauffman net. Because the Kauffman net has many fewer relevant 
modules than the lattice system, even for large N, one will be able to 
observe realizations which always relax into a time-independent be- 



havior, a cycle of length one. According to Bastolla and Parisi [1998b] 
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the average number of cycles observed in a given realization will be 
independent of N. Realizations with long cycles will be exponentially 
unlikely. However, when one gets a long cycle, one can expect to have 
very many of them. The frozen phase is one in which almost all the 
variables end up with fixed values. That is, for a given realization, 
after a large number of steps most variables approach a value which 
is independent of the initial state of the system. Thus the system be- 
haves, for most of the variables, as if it were highly frictional with the 
result that each variable comes to a stop its own 'best' value. This 
is the reason that the overlap goes to one as time goes to infinity in 
this phase. 



The majority of our knowledge of the Kauffman-net behavior of the 
frozen phase comes from two calculations. As discussed above, the Kauff- 



man net with K = 1 was solved exactly by Flyvbjerg and Kjaer (Flyvbjerg 
and Kjaer [1988 1 ) . For the case K = 2, Lynch [1993h| ] has proven a group 



of theorems which apply on the frozen side of the K = 2 system. To en- 
sure the system was frozen, Lynch looked at the case in which there is 
unequal weight to the four classes of K = 2 functions, and demanded that 
the constant functions have a larger weight than the non-canalizing ones. 
This ensured that the Hamming distance, as calculated by an approach like 



that in equation (2.4), would decay exponentially. He then proved that the 



system was in a frozen phase by showing: 

1. Almost all elements in the system were in their stable state. 

2. Almost all gates were weak, that is, changing their value does not 
affect the cycle that is entered. 

3. The length of the transient, i.e. the number of steps before the system 
enters its cycle, is bounded below by a constant times In TV. 



4. There is also a rather strong bound on the cycle length. The bound 
includes the statement that the median cycle-length is bounded by a 
number which is independent of N. 



3.4.2 Chaotic Phases 

In contrast to the frozen phase, the chaotic phase is one in which a nonzero 
fraction of the variables remain sensitive to initial conditions. In fact, most 
variables keep changing their values throughout the time development of 
the system. In the chaotic phase the average length of limit cycles and of 
the transient which occurs before the entry of the limit cycle both grow 
exponentially with N. Because of this explosive growth as the system is 
made larger, numerical investigations of orbits in the chaotic phase are 



limited to small system sizes (see, e.g., Bastolla and Parisi [1998b]) 
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As discussed above, the case when each input is coupled to all the others 
can be solved exactly (Dcrrida and Pomeau [1986 1; pcrrida and Flyvb 



jerg [1987b||), with the typical cycle length growing with system size as 



2 N / 2 . When K is finite but greater than the critical value K c , both the 
Kauffman net and the lattice models have typical cycle lengths that grow 
exponentially with N. This behavior reflects a complex network of depen- 
dency loops in this regime. Though it is plausible that the details of the 
interconnected loop structures could be different for lattice models and for 
the Kauffman net, they do not appear to lead to marked differences in the 
behavior of the two types of models in the chaotic regime. 



3.4.3 Critical Behavior 



We have already seen that there is a close connection between the phase 
transition of the N-K model and the standard percolation transition of sta- 
tistical physics. According to Obukhov and Stauffcr [ 1 9 89| ] , random Boolean 
networks on nearest-neighbor d-dimcnsional lattices belong to the 'univer- 
sality class' of directed percolation with quenched disorder in d+1 dimen- 
sions. However, this result is an argument, not a theorem. The sharing of 
the universality class means that the N-K models share many of the de- 
tailed properties of correlation and ordering with the directed percolation 
models of standard statistical mechanics. More particularly, the annealed 
N-K model based upon nearest neighbor interactions is equivalent to a di- 
rected percolation problem ( Derrida and Stauffcr [1986| ). Moreover, when 
the interaction range is infinite (roughly corresponding to infinite dimen- 
sions or a simplectic geometry), the quenched (or usual) N-K model is 



equivalent to the corresponding directed percolation problem (Dcrrida and 
Pomeau [1986| ; pcrrida and Weisbuch [1986[ ) 



Many studies have been carried out to determine the statistical properties 
of cycles in the critical phase. The main results can be summarized by 
saying that in the critical phase, both the typical cycle lengths as well as 
the mean number of different attractors grow algebraically with N . More 
explicitly, it was believed that these quantit i es were proportional to y/N 
see for instance, Bastolla and Parisi [1998b| ; |Flyvbjerg [1989| ; |Kauffman| 



[1990 



196S 



1995 



199S ] ) . This was one of the most attractive results of the 



N-K model in that it matched one of the p owe r-law behaviors exhibited 
by living organisms, as we will see in section 5.1. However, S. Bilke and F. 
Sjunnesson have recently brought into question the y/~N dependence of the 
mean number of different attractors, arguing that this result comes as a 



consequence of a biased under sampling of the whole state space (Bilke and 
Sjunnesson [200l| ). 



Because the size of the state space grows exponentially with the number 
of elements N, it is not feasible to enumerate completely all the different 
cycles in a given network realization unless N is quite small. What people 
usually do is to probe the state space with a small fraction of its elements 
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FIGURE 3.5. Mean number of different limit cycles as a function of the network 
size for critical Kauffman networks with connectivities K — 2, 3, 4. The solid 
curves are the best linear fits to the numerical data, whereas the dashed curve is 
a v/iV fit for the case K — 2. Taken from 



Bilke and Sjunnesson [2001] 



(which are assumed to be representative), and then to infer the statistical 
properties of the model from this sampling. However, this method, which 
has yielded the y/N dependence referred to above, has the problem that it 
is possible to miss a small but still important fraction of cycles. 

In their work, Bilke and Sjunnesson (ibid.) use a different approach. They 
present a decimation method which eliminates the stable elements of the 
network^], retaining only the relevant ones. Since for any realization of the 
network the number of different cycles as well as their lengths depend only 
on the relevant elements, all the statistical information related to limit 
cycles is preserved through this decimation procedure. The advantage in 
eliminating the stable elements is that the number of variables is drastically 
reduced. As a consequence, it is much easier to perform a full enumeration 
study of the number of different cycles for each realization of the network. 
Through this approach, Bilke and Sjunnesson found that, in the critical 
phase of the Kauffman net, the mean number of different attractors (iV c ) 



8 A stable element is one which evolve to the same fixed value independently of the 
initial state. 
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grows linearly with N instead of as v N. Fig. |3.5| (taken from Bilke and 



Sjunnesson [2001 1) shows the mean number of different cycles (N c ) as a 
function of the system size N. As can be seen, the linear dependence (N c ) ~ 
N fits the numerical data much better than the y/~N behavior. 

In addition to the exponential growth of the state space as N increases, 
the number of network realizations grows supcrexponentially with K . The 
evidence presented by Bilke and Sjunnesson raises the possibility that some 
of the results obtained so far related to cycle lengths and basins of attraction 
could also include a systematic bias due to undersampling. More work is 
needed to address these questions. 



4 Reversible Models 

Kauffman nets are generic models for N elements coupling K different 
variables. In this section we specialize this generic model to include an 
interesting symmetry property: time reversal symmetry. The discrete sym- 
metry and quenched randomness together cause some new properties to 
emerge. 

The standard Kauffman net that we have considered so far, whose dy- 



namics is given by equation (1.1), is dissipative because multiple different 
states can map into one, so that information is lost. In this section we 
will refer to this system as the dissipative Kauffman net or simply as the 
dissipative N-K model. Not all systems are dissipative; many of the sys- 
tems considered in Hamiltonian mechanics are reversible. Reversible sys- 
tems have the property that some transformation of the coordinates (for 
example changing the sign of all velocities) makes the system retrace per- 
fectly its former path. This section discusses some aspects of the behavior 
of discrete reversible maps. In the time-reversible Boolean network studied 



Coppersmith, Kadanoff, and Zhang [2001a 0, the state of the system at 



time t + 1 is governed by the equation 

o*(t+ 1) = ■ ■ ■ © <Ti{t - 1), (4.1) 

where the © denotes addition modulo 2. Each time-reversible network real- 
ization has a corresponding dissipative realization with the same functions 
and connections. In the dissipative N-K model, the state St+i is com- 
pletely determined by the previous state E t . Nevertheless, from equation 



(4.1) we see that in the reversible model both £t_i and E t are required 
to calculate £4+1. Thus, in the reversible model the state of the system at 
time t is represented as 



S t = 



and the state space has now 2 2N points. The behavior of the reversible 
model is in some ways closely analogous and in other ways quite different 
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than that of the dissipative model. Both models exhibit a phase transition 
at K c , a critical value of K below which the system "breaks up" into 
disconnected sections, and above which there is a percolating cluster of 
connections. The value of K c is slightly lower for the reversible model (K c 
1.6) than for the corresponding dissipative model (K c = 2), and some 
details of the transition are different. But for both models the observed 
behavior is consistent with a percolation picture. For K < K c , changing 
one element of a system leads to changes in only a finite number of other 
elements, whereas for K > K c , changing one element causes a cascade of 
influence that spreads to a nonzero fraction of all the other elements. The 
typical cycle length grows slower than linearly with the size of the system 
when K < K Cl and exponentially with system size when K > K c . 

One big difference between the reversible and dissipative models is that 
they have vastly different numbers of attractors (or, in the reversible case, 
limit cycles) . Large differences in the behavior are not entirely unexpected 
because in the dissipative model many different state space points end up at 
the same attractor, whereas in the reversible model every state space point 
is on exactly one limit cycle. For example, when K = (all input functions 
either 1 or 0), the usual Kauffman net has only one attractor, while the 
number of limit cycles of a reversible Kauffman net is proportional to 2 2N . 
The reversible result can be understood by noting that for either input 
function, each element is in one of four different cycles, depending on its 
initial conditions, and that when K = the elements are all independent. In 
the other limiting case K ~ N, as discussed above, the number of attractors 



of a dissipative Kauffman net is proportional to N (see equation (3.9)). In 
contrast, in the reversible model with K = N, and indeed throughout the 
regime K > K c , the number of limit cycles grows as 2 N . 

The number of attractors in the reversible model with K = N can be 
understood by studying the mechanisms that lead to orbit closure. One 
way to close an orbit is to repeat two successive configurations, so that 
Yit = So and Ey + i = Ei. Using the approximation that each successive 
E is chosen randomly from the 2^ possibilities yields a probability of an 
orbit closure at a given time of order (2~ N ) 2 = 2~ 2N . However, there is 
another mechanism for orbit closure that leads to much shorter orbits and 
thus totally dominates the behavior. If at some time r one has E T = E T _i, 
then time-reversibility implies that E T+ „ = E T _i_„ for all n. Similarly, if 
E r+ i = E T _i, then E T+rl = E T _„ for all n. Because the orbits reverse at 
these points, we call them "mirrors." Once two mirrors have occurred, then 
the orbit must close in a time that is less than the time it has taken to hit 
the second mirror. Again assuming that each successive E is chosen ran- 
domly from the 2 N possibilities, one finds that a mirror occurs at a given 
time with probability proportional to 2~ N , so the expected number of steps 
needed to hit two mirrors is of order 2^. Hence, typical orbit lengths are 
of order 2 N . Since there are 2 2N points in the state space altogether, the 
number of limit cycles is proportional to 2 N . In the reversible Kauffman 
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net, when K is finite but greater than the critical value K c , the distribu- 
tion of cycle lengths can be extremely broad. For example, when N = 18 
and K = 2, the median cycle length I is approximately / m 140, and yet 
the probability of observing an orbit of length 10 10 is greater than 10~ 4 . 
This huge variability arises because of a nontrivial interplay between the 
discrete symmetry and the quenched randomness in the system. The occa- 
sional extremely long orbit arises because some combinations of coupling 
realizations and initial conditions are such that mirrors cannot occur at all. 
If the mirrors are not available to close the orbits, then the system must 



wait until two successive configurations happen to repeat. In (Coppersmith 



Kadanoff, and Zhang [2001a |) it is shown that in a system with finite K in 
the true thermodynamic limit, almost all realizations and initial conditions 
yield no mirrors, and typical orbit lengths grow as 2 2N ( 1 ~ e \ where e is of 
order l/(2 2 ). However, the crossover to the limiting behavior occurs only 
when N ~ 2 2 , so that even for moderate K, this behavior is not accessible 
numerically (for example, when K = 3, one requires N ~ 256, a value at 
which 2 N ~ 10 77 and 2 2N ~ 10 154 ). 

The significance of this enormous variability in the cycle lengths is not 
clear. One way to interpret these results is to conclude that characterizing 
cycles is not the right way to study the model. It would be interesting to in- 
vestigate whether enormous variability arises in other random systems with 
discrete symmetries, and to determine whether there arc possible experi- 
mental consequences. One possible starting point for comparison is with 
properties of random mag nets, for which regions of atypical couplings lead 
to Griffiths singularities ( Griffiths [1969 1 ) , which have stronger effe cts on 
dynamic than static properties ( Randcria, Scthna, and Palmer [1985 1 ) . One 
must explore whether the analogy is appropriate (this is not obvious be- 
cause Kauffman networks are not lattice-based), and if so, whether the 
results for the reversible Kauffman model indicate that dynamics far from 
equilibrium can be even more strongly affected by atypical coupling real- 
izations than are the properties in spin models that have been studied to 
date. 
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5 Beyond Cycles 

5.1 Non-robustness of Cycles 

A cycle forms when the system returns to a configuration which is exactly 
the same as one it has previously visited. A demand for an exact return 
should be viewed not as a single constraint but in fact as N constraints 
upon the configuration. Such a strict demand makes the properties of cycles 
quite special and delicate. For this reason, a study of cycles is probably not 
what one would want for understanding the possible physical or biological 
consequences of the models like the N-K model. The statement just made 
flies in the face of a very large body of work, some of which we have just 
described. We should, for this reason, argue for this statement with some 
care. 

Why don't we believe in cycles? 

1. The characterizes of cycles are neither intensive nor extensive vari- 
ables. 

2. With exponentially short cycles in localized regions one gets power 
laws overall. 

3. The critical situation has very many, short cycles which are not ob- 
served when one starts from randomly picked starting points. 

4. Attractor basins are complex in character, being at best multi-fractal 
for large chaotic systems. 

5. In a large system you must wait so long to see a cycle that it cannot 
be really important. 

6. In a large chaotic system, changing one rule changes the cycles quite 
a bit. 

7. In both glasses and biological systems one wants to characterize the 
system by things which are very robust. 

We do believe that generic properties of cycles are important though, in 
that they characterize general aspects of dynamical systems. Nonetheless, 
the huge fluctuations throughout realizations in such important quantities 
as cycle lengths and number of different attractors, calls for other types of 
characterizations. Real networks, whether they are genetic or neural or of 
any other kind, arc always subjected to external perturbations. The robust- 
ness in the dynamics of the network can not rely on quantities which change 
dramatically with perturbations. Hence, it is important to characterize the 
dynamical properties of the network in the presence of noise, trying to find 
out which kind of quantities are preserved under the influence of noise, and 
which ones are not. 
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5.2 Characterization via noise 

The addition of noise to a map provides a possibility for generating addi- 
tional information about the behavior of the models. Noise naturally blurs 
out the sharpness of behavior, making for a more "fuzzy" characterization. 
Noise is, then, a natural way to get away from the difficulty posed by the 
overly-precise characterization provided by the cycles. 

Unfortunately, most previous work on the N-K model does not include 



noise. We do wish to poin t to two papers ( Miranda and Parga [1989 ] 



Golinclli and Dcrrida [198S ]) in which noise has been used to probe N-K 
behavior. For reasons which will become more evident later on, we describe 
these papers respectively as a "crossing paper" and a "convergence paper" . 
These papers both break the precision of the dynamical rule of equation 
by saying that the rule is broken with a probability r 



<7<(*+l) = 



fi i u h (i) (*)>•••! a m i (i) (*)) with probability 1 
"■/« {°3i (») (*)>•••) a i Ki (i) (*)) witn probability 



(5.1) 



which can also be written in the alternative form: 

{fi ( a h (i) (*) ! • • • i a m i W (*) ) witn probability 1 - 2r 
1 with probability r 

with probability r. 

(5.2) 

This equation can be described as providing probabilities r for the two 



possible values of the outcome, independently of the value of /j. In Golinclli 



and Derrida [1989 ] the rules are described in terms of a temperature T, 



related to r by 

1 - tanh(l/T) 



(5.3) 



Both groups examine the development of two or more different initial 
configurations using the same realizations. They also apply exactly the 
same rules and the same probabilistic choices to the different configura- 
tions. So far, both papers are essentially similar. There arc two kinds of 
differences, the first being the choice of measurement, and the second being 



the way they apply equations ( |5.l| ) and (p.2j). 

The convergence paper starts with two or more randomly chosen initial 
configurations, Ej, £§,..., E™, and calculates the resulting trajectories 
step by step: 

v 1 — > v 1 — > v 1 — > v 1 

n ^ n l ^ n 2 ^ ■ ■ ■ n T 
v^2 v^2 

^0 — * ^1 ~~ * ^2 ~~ * • ■ ' 



srrm y-mi vim 

Zjq — * Zjj — > Zj 2 — + . . . Zj r 
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At each step, and for each i, a choice is made among the three branches 
of equation (5.2), and that choice is equally applied to the m configura- 
tions.^ The calculation continues until two of the m configurations become 
identical (say for example E* = E^). The time needed to achieve the con- 
vergence is noted. We will denote this time by r m , stressing the fact that 
m configurations are being analyzed. 

In some ways, the convergence calculation is mor e complicated than the 
crossing one. The noise as defined by equation (5.2) tends to produce con- 
vergence because it is applied equally to all trajectories and because makes 
the values of the elements to be equal in all trajectories. In the limit of 
infinite N, the system shows three phases: the low noise phase in which 
almost always trajectories will not converge, the high noise phase in which 
trajectories will always converge, and a separating critical phase. These are 
respectively described as low temperature, high temperature, and critical 
phases. 

The authors examine K = 4 and large- N, picking the noise-levels r = 
0.15 (low temperature phase), r = 0.25 (critical phase), and r = 0.35 (high 
temperature phase). Without noise, the N-K model would show chaotic 
behavior, and very long cycles. In the presence of noise the average of 
lnr, with r being the convergence time, diverges as N goes to infinity. The 
behavior of the convergence time for two orbits, T2, can be calculated in the 
annealed sequential update case, and the results compared with simulations 
for both annealed and quenched system. In both simulation and theory, the 
three phases are characterized by having different forms of divergence of 
(I11T2) with N (see Fig. 5.1). At low temperatures the divergence is linear: 
(I11T2) ~ N. In the critical phase (I11T2) diverges linearly with InN as 
(lnT2) = 0.5 In N. The weakest divergence occurs in the high temperature 
phase in which (I11T2) varies as lnlniV. These A^-depcndcncics describe 
the variation of the number of elements forming the barrier to having two 
configurations merge into one another. The noise causes these elements 
to be identical and then the merge occurs. Thus, for example, in the low 
temperature phase, we must bring to equality a finite fraction of all the 
elements in the system in order to have the convergence. 

In some loose sense, these numbers measure the size of the barriers which 
hold together the attractors for this system. But it is hard to know what 
the attractors themselves might be. Because r is relatively large, these 
presumed attractors are probably not the cycles of the original system. In 
fact, the behavior of T2 is much the same for the annealed system (which 
has no cycles) as for the quenched system, which does have cycles at r = 0. 
Nobody has yet investigated the limiting case as r goes to zero. It might 
be most interesting to look at this limit, particularly in association with a 



9 The two-branch versus three-branch methods become incquivalent when they are 
applied to several configurations at once. 
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FIGURE 5,1. Plots of (lnr m ) versus N, taken from reference Golinelli and Per 



rida [1989 1. In all the cases the connectivity of the network was K = 4. The three 
graphs correspond to three values of r: (a) Low temperature phase, r = 0.15; (b) 
Critical phase, r — 0.25; (c) High temperature phase, r = 0.35. 



limit which keeps the system critical (say K goes to 2). 

Let us consider the case in which more than two trajectories are analyzed. 
Start with m trajectories and let r m measure the first time when any two 
of these have converged. If there are only M large or important basins, one 
might well expect (lnr m ) <C (In tm) for m > M . Instead one observes that 

(hiT m ) = (lnr 2 ) - In , 

in both theory (in the annealed approximation) and simulation. This form 
indicates an indefinitely large number of attractors, all with basins of com- 
parable size. 

The crossing paper considers two different configurations and Sg, 
which can belong either to different attraction basins or to the same basin 
of attraction, and then iterates forward, noting all configurations {Sj, Y\, 
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FIGURE 5.2. (a) Log-log plot of the crossing time r as a function of the level of 
noise r. The solid curve is th, the crossing time starting from different attractors, 
whereas the dashed curve corresponds to t s , the crossing time starting from 
the same attractor. (b) The index A(r) describing the fractal dimension of the 
attractor plotted as a function of r. The solid curve is the result of the numerical 
simulation, obtaine d as A(r) — 21nr/lnf2, and the dashed line is the theoretical 
prediction given in (5.7). Note that these curves become identical when r — > 0.5. 
Both graphs (a) and (b), calculated for a situation with K = 2 and N = 20, were 
taken from reference |Qu, Kadanoff, and Aldana [2002], 



• ■ • } an d {S 2 ,, ^li ^ii ■ ■ ■ } they produce. In each step of iteration, and 
for e ach i- value, a choice is made between the two branches of equation 
(5.1), and that choice is applied to both configurations. This continues 
until the time r in which one of the two trajectories attains a configuration 
previously entered by the other one (for example, if £ 2 is equal to one of 
the configurations {£j, £}, . . . , £*}). The calculation is terminated at this 
crossing event. The measured quantity is the "time" t needed to produce 
the crossing. 

In the absence of noise, if two initial configurations belong to different 
basins of attraction, the time for the two subsequent trajectories to cross 
is infinite. In the presence of noise, there is a chance for each trajectory to 
"jump out" of its basin of attraction, exploring a bigger part of the state 
space. The two trajectories will have a number of opportunities equal to 
r 2 , to cross one another before r steps have elapsed. If the size of the space 
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being explored by the trajectories is f2(r), then the typical time for the 
crossing will be 

T«ft(r) 1/2 . (5.4) 

Miranda and Parga simulated the system and measured r as a function of 
r and N for the critical Kauffman net in which K = 2. Their result for 
large N may be summarized as 

n(r) = n A(r \ (5.5) 

where = 2 N is the volume of the state space in the system. 

In the work by Miranda and Parga the "fractal" exponent A(r) was 
not estimated accurately. This work was recently extended by X. Qu et 
al., who consider larger values of noise and different connectivities of the 
network ( Qu, Kadanoff, and Aldana [2002f |). The authors analyze two cases 



to compute the crossing time, when the two initial configurations belong 
to the same basin of attraction, and when they belong to different basins. 



We will denote these two crossing times by t s and n respectively. Fig. 5.2a 
shows the average crossing times Td and r s as funcions of r, for a net with 
N = 20 and K — 2. As can be seen, when r is close to its maximum value 
0.5, both times are practically the same. In fact, Qu et al. have shown that 
for large values of r both t s and behave as 



The above expression agrees with equations (5.4) and (pJ5J) by identifying 



In [1 + (1 - 2r) 2 /2 K ] 
Mr) = 1 - - 1 \ n2 1 ■ (5-7) 

In contrast, when r is close to the behavior of Td and t s differ substan- 
tially. For r — > the divergence of Td is given simply by 

Td^d/r + TQ (5.8) 

where C\ only depends on K and N, and To is the value of t s at r = 0. 
The complete analytical cxpresion of A(r), valid in the whole interval 



[0, 1/2] is not known yet. Fig. 5.2b shows a plot of A(r) obtained by nu- 
merical simulations. It is interesting to note that the attractor has a fractal 
volume which depends upon r. Once again, one is frustrated because one 
does not know what the attractor might be. It is once again probably not 
anything directly related to a cycle, since starting points in the same cycle 
or in different cycles both give the same r- values for the higher values of 
N and r. Here too one might guess that studies with smaller values of r 
might shed light on the N-K model attractors. 
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6 Applications 

6.1 Genetic networks and cell differentiation 

N-K models have been widely used in the modeling of gene regulation 
and control ( Somogyi and Snicgoski [I996fl; Dc Sales, Martins, and Stariolo 



|[1997| ; [Kauffman [1993| ; [Huang and Ingbcr [2000| ). A very remarkable 



characteristic of multicellular organisms is that all the different cells of 
which they are made have the same genetic information. What makes the 
difference between the different cell types are the genes which are expressed 
in every cell at every moment. In a given cell type, some particular genes 
are turned off and others are turned on. So, in a liver cell, only the "liver 
genes" arc being expressed while all the other genes are turned off, whereas 
in a neuron the "neuron genes" are the only ones which are expressed. 

The physical and chemical mechanisms by which the cell determines 
which genes are to be expressed and which are not are not yet fully un- 
derstood. There is evidence that gene regulation and control can occur at 
every stage during the metabolic pathways leading up from the genetic 
information contained in the DNA, to the translation of this information 
into proteins. Nevertheless, most of the gene regulation and control seems 
to occur at the level of transcription of the genetic information. At this 
level, one gene of DNA is transcribed into a molecule of messenger ARN 
(mARN) only if the conditions for this transcription are present. In the 
most simple model (applicable to bacteria), for the transcription of one 
gene into a molecule of mARN, it is necessary a protein, called activator, 
which attaches to the beginning of the gene indicating that this gene is 



ready to be transcribed (sec Fig. 6.1a). On the other hand, there also ex 



ist repressor proteins which, when attached to the beginning of the gene, 
inhibit its transcription, turning the gene off (see Fig. |6~l|b). 

In eucaryotic cells the situation is more complicated in that many acti- 
vator or repressor proteins might be needed to activate or to inhibit the 
expression of a single gene. For example, it is known that the human /3- 
globine gene (expressed in red blood cells) is regulated by more than 20 
different proteins. Some of these proteins may function as both repressors 
or activators, depending on how they are assembled. 

The activator or repressor proteins of a given gene are themselves speci- 
fied by some other genes, whose expression, is in turn controlled by other 
proteins codified by other genes and so on. Genes interact each other 
through the proteins they specify: the product protein of one activated 
gene can influence the activation or deactivation of other genes. Similarly, 
the absence of the product protein of a deactivated gene can influence the 
activation or deactivation of several other genes. In some particular cases, 
it is known which gene controls which other one, but in most of the cases 
the interactions between genes in a given cell type are completely unknown. 
A picture emerges in which genes are mingled together forming a network 
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FIGURE 6.1. Schematic representation of gene regulatory proteins, (a) An ac- 
tivator protein attaches to the gene at a specific binding site, activating the 
polymerase II which then can transcript the gene into a molecule of mRNA. (b) 
When a repressor protein is attached to the gene, the polymerase II is blocked 
up and therefore unable to transcript the information contained in the gene. 



of interacting elements which are connected in a very complicated way. 

The N-K model was first suggested by Kauffman as a way of model- 
ing the dynamics and evolution of this complicated network of interacting 



genes ( Kauffman [196S ]). Within this model, two genes are linked if the 
product protein of one gene influences the expression of the other one. In 
real cells, genes are not randomly linked. Nevertheless, the web of linkages 
and connectivities among the genes in a living organism is so complicated 
(and mostly unknown), that to model these genes as randomly linked seems 
to be a reasonable approximation. By this means, the results coming out 
of the model are not restricted to a particular set of linkages and connec- 
tivities, but reflect instead the generic properties of this class of networks. 

The state of the cell at every moment is determined then by the state 
of its genes a±, . . . , (Jn ■ Before proceeding further, it is important to recall 
the main assumptions that are usually made to describe the dynamics of 
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genetic networks by means of the N-K model: 



Every gene o~i can be only in one of two states, turned on (1) or 
turned off (0). 



• Every gene is randomly linked with exactly K other genes. 

• The evolution rule fi associated to gene Oi is a weighted function 
which acquires the value 1 with probability p ant the value with 
probability 1 — p. 

• The updating of the network is synchronous. 

As we have seen, under these assumptions the phase space breaks up 
into different cycles or attractors whose properties in the frozen phase and 
the chaotic phase are substantially different. According to Kauffman's in- 
terpretation, each attractor represents a cell type or a cell fate, whereas a 
single state of the system represents just a temporal state which the cell is 
passing through. From this point of view, cell reproduction processes start 
with an initial configuration of genes, which eventually evolves towards its 
corresponding attractor; the attractor determines a particular cell type or 
cell fate. 

For this mechanism of cell differentiation to be meaningful, the length 
of the attractors must not be too long, for otherwise the cell would never 
reach its stable cycle. In the chaotic phase the length of the cycles grows 
exponentially with the system's size (I ~ 2 aN , where a is of order 1). 
Therefore, the system has to go through very many states before reaching 
a stable set of configurations. In addition, systems in the chaotic phase are 
very sensitive to perturbations (mutations), partially because the number 
of relevant elements is comparable to the size of the system. The above 
prompts the thought that genetic networks of living organisms are not in 
the chaotic phase. 

On the other hand, in the frozen phase the cycles are much shorter 
than in the chaotic phase: an initial set of genes swiftly reaches its stable 
configurations. The fraction of relevant elements in the frozen phase is 
close to zero. As a consequence, the system is extremely resistant to point 
mutations |^] or to damage in one or more of the evolution rules. But in 
order to evolve, genetic networks of living organisms should allow some 
degree of sensitivity to mutations, which rules out the frozen phase as a 
physical state which living organisms could be in. 

Kauffman suggested that gene networks of living organisms operate at 



the edge of chaos (Kauffman [1993 1 ) , meaning that the parameters have 
been adjusted through evolution so that these networks are at or near the 
critical phase. In the critical phase, both the number of different attractors 



10 A point mutation is a change in the value of one gene. 
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as well as their lengths are proportional to a power of N, so the cell can 
reach very quickly its stable configurations. Also, the fraction of relevant 
nodes in critical networks, even though is small, is not zero, which means 
that these kinds of networks present some degree of sensitivity to changes in 
the initial conditions. In other words, critical networks exhibit homeostatic 
stability, a term which we will come to in the next section. 

A very remarkable observation supporting the idea of life at the edge of 
chaos consists in the fact that the number of different cell types in an organ- 
ism is roughly proportional to the square root of its DNA content. Further- 
more, the mitotic cycle period, which can be considered as a measure of the 
time required for a cell to reproduce, seems also to be proportional to the 
square root of the cell's DNA content ( Kauffman [1993| ). Thus, random 
networks in the critical phase seem to satisfy the requirements of order, 
evolvability and stability found in living organisms. 

Even though this idea is very attractive, there are some problems yet to 
be solved. As we have mentioned, for unbiased evolution rules /s, the criti- 
cal phase is characterized by the low connectivity K = 2. This implies that 
genetic networks of real organisms are restricted to have very low connec- 
tivities in order to be at "the edge of chaos" . As soon as the connectivity 
grows, the system becomes more and more chaotic. But it is well known 
that the connectivity in real genetic networks is rather high. For example, 
the expression of the even-skipped gene in Drosophila is controlled by more 
than 20 regulatory proteins, also the Human /3-globine gene we have re- 
ferred to before ( Alberts, Bray, Lewis, Raff, Roberts, and Watson [1-H)j[ ) . 
In eucaryotes it is common to find that one single gene is regulated by a 
bunch of proteins acting in association. On the other hand, sometimes when 
a single signaling receptor protein is activated, it can influence directly the 
activation of a very large array of genes. Let us consider for instance the ac- 
tivation of the platelet-derived growth fa ctor (3 receptor (PDGFR/3), which 



induces the expre ssion of over 60 genes ( Fambrough, Mcclure, Kazlauskas 



and Lander [199S Q. These examples, among many others, suggest that the 
connectivity in real genetic networks is not low, but on the contrary, it is 
very high. Nonetheless, cells do not seem to operate in the chaotic phase. 

There are two ways to increase the connectivity in the N-K model with- 
out going out of the critical or ordered phase: (a) by the use of weighted 
evolution functions /j, or (b) by the use of canalizing functions. When 
evolution functions fa are weighted with the probability parameter p, the 



critical line is given by equation (2.8). Fig. 2.4 shows the graph of K c as a 



function of p. Even if we suppose that genetic networks can be either in the 
frozen phase or in the critical phase, they would be restricted to remain 
within the shaded area of the figure. 

On the other hand, it is known that the fraction of canalizing functions 
becomes very small as K increases (Kauffman [1984]). An upper bound for 
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this fraction is 

AK 



2 21 



(6.1) 



which tends to zero as K — > oo. Consequently, canalizing functions are 
extremely rare when K is large^j. If the apparent order seen in living cells 
relies on either weighted or canalizing functions, it has still to be solved 
what kind of mechanisms drove, through evolution, the genetic networks 
towards the generalized use of such type of functions. 

Another problem lies in the assumption that genes can be in only two 
states. It is true that a given gene is expressed or is not. But the prod- 
uct protein of a gene can participate in a variety of metabolic functions, 
producing even opposite effects depending on the physical and chemical 
context in which it acts. Such behavior can be modeled by assuming that 
every gene can acquire more than two states. But in such a case the con- 
nectivity of the network must be even smaller to keep the system within 
the ordered phase. For if we assume that every gene can be, on average, in 
one of to possible states, and if every one of these states is activated with 
the same probability, the critical connectivity K c is then given by 

K c =^-. (6.2) 

TO — 1 

The critical connectivity decreases monotonically when to > 2, approaching 
1 as to — > oo. The moral is that for this kind of multi-state networks to 
be in the ordered phase, the connectivity has to be very small, contrary to 
what is observed in real genetic networks. 



Partially to overcome these difficulties, Dc Sales, Martins, and Stariolo 



[1997 1 have proposed a model which is slightly different from the original 
N-K model. In their model, the dynamics of the system is governed by the 
equation 



*i(t + 1) = Sign | J im {t) + J «iW ii(i)(*)| ' ( 6 ' 3 ) 

where Jij t u) is the coupling constant representing the regulatory action of 
the input of gene i, (I = 1, 2, . . . , K — 1). There is also an autogenic 
regulation, given by J a. The set of coupling constants Jy represents the 
very complicated set of biochemical interactions between genes. Since these 
interactions are mostly unknown, the authors assign the coupling constants 
in a random way, according to the following criteria: 

• The product protein of a given gene can activate, inhibit, or not affect 
at all the transcription of another gene. 



11 Even for K as small as K = 10, equation ( |6.l| ) gives a fraction of canalizing functions 
of the order 10~ 153 . 
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FIGURE 6.2. Phase space for the N-K model with the dynamics given by 



equation (3.4). As can be seen, under this dynamics the system exhibits three 



different "phases": frozen, marginal and chaotic (taken from 
and Stariolo [1997| ). 



De Sales, Martins 



• Interactions between genes are not necessarily symmetrical, namely, 
Jij 7^ Jji in general. 

• Autogenic or self-regulation is allowed. 

The coupling constants are then assigned according to the probability 
function given by: 

P(Jij) = l —^ L Wa - 1) + S(Jij - 1)] +PiS(J ij ). (6.4) 

In this way, the couplings Jij can be activating (+1) or inhibitory (—1), each 
with probability (1— pi)/2, or neutral (0) with probability pi . Furthermore, 
the linkages for every gene are chosen either at random among the whole 
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set of genes, with probability P2, or only among the nearest neighbors, with 
probability 1 — By varying the parameter -pi one can go from lattices 
where local interactions are the most important, to random nets where all 
ranges of interaction are present. 

By analyzing the Hamming distance between two initial configurations, 
the authors show that there are three phases in the dynamics of the system: 
a frozen phase, a marginal phase and a chaotic phase, as shown in Fig. |6.2| . 
The frozen and chaotic phases are as in the traditional N-K model. The 
marginal phase is characterized by the fact that all the attractors are stable 
since a small change in the initial configuration neither vanishes nor grows. 
Also, in this phase both the number of different attractors and their lengths 
grow as a power of N. The above results were obtained for networks with 
connectivity K = 9 and N = 625. It is probable that the marginal phase is 
only a finite-size effect, vanishing for an infinite system. More work has to 
be done to explore the whole space of parameters of this model, but these 
preliminary results show that it is possible to obtain ordered and marginal 
behaviors even in the case of high connectivities by changing the type of 
interaction between the genes. 



6.2 Evolution 



In a traditional framework of the theory of evolution ( Simon [1969( |; Dawkins 



[19891 , |l986l ]), chan ges in the phenotype of organisms are accumulated grad 



ually, yielding a gradual increase in complexity of form and function. Sim- 
ple structures slowly assemble together to form more complex structures, 
which in turn assemble to build up even more elaborate systems, and so 
on. At every stage in the formation of a complex system (organism) out 
of simpler elements, many sub-systems are created, which represent tem- 
poral stable states along the way in th e construction o f the whole system. 
As T. D. Zawidzki has pointed out in Zawidzki [1998|1 , from this point of 



view the problem of evolution of complex systems translates into a search 
problem. For the construction of a complex system, evolution searches in 
the space of the possible configurations leading to the whole system. If 
no stable intermediate configurations were found in this search, the con- 
struction of the whole system would be extremely improbable, since the 
number of possible configurations in which the parts can be arranged in- 
creases exponentially with the number of different parts. Evolution does 
not search in one step for the "final" configuration of a complex system. 
Instead, the search is carried out gradually, finding first intermediate stable 
configurations of sub-systems that are then assembled, giving rise to the 
whole complex system. Every sub-system solves a particular problem (or 
set of problems) , facilitating the construction and functioning of the whole 
organism. Furthermore, once an evolutionary problem has been solved by 
primitive organisms, more complex species which evolve from them still 
continue solving this problem in the same way. New species are faced with 
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new problems, but still preserving the old solutions to the old problems (the 
genetic code, for example, was "invented" only once). According to this the- 
sis, genomes of complex organisms are made up of functional modules of 
genes, each module encoding the solution of a given biological problem 
encountered by the species at some point through evolution. 

The study of Boolean networks has suggested new mechanisms for evolu- 
tionary processes. Kauffman has stressed that many evolutionary changes 
involve reorganizing of the same genetic material rather than making it 
more complex flKauffman [1993 , 1995 ]). In Kauffman's approach, genes are 



organized in complex genetic networks provided with some given dynamics. 
The "search" of evolution consists in searching for the most stable orga- 
nization of genes, which in terms of the N-K model means stable cycles. 
From this point of view, there is not an increasing complexification of the 
evolving system, but a inherent complex organizational dynamics which 
settles down in a finite number of stable attractors. The "role" of evolution 
then is to look for the more stable attractors which the system can fall into. 

A Boolean network made up of N genes has 2 N states, but the system 
organizes itself into a much smaller number of attractors. Depending on 
the parameter values being used (K and p) , these attractors are stable or 
they are not. Stability is defined according to the response of the network 
to perturbations, which can be of three different kinds: 

• changes in the states of a few genes by flipping the value of some 
randomly chosen ones; 

• permanent changes in the linkages of some genes; 

• permanent changes in the values of the evolution functions fi associ- 
ated with some genes. 

As we have already mentioned, only networks in the critical phase have 
the stability required to constitute evolvable systems, in that these net- 
works are able to recover to most of the mutations described above. In the 
chaotic phase, the attractors are extremely unstable since any kind of per- 
turbation would shift the system to another attractor. In the frozen phase, 
even though the majority of the genes are motionless, small changes in the 
linkages or in the evolution functions would make the system jump to a 
very different attractor if these damages are carried out on the relevant 
elements of the net. But networks operating in the critical phase show a 
very high homeostatic stability, which means that after some perturbation, 
such networks are very likely to fall again in the same attractor. Conse- 
quently, in this approach, evolution is also interpreted as a sort of search 
for gene networks possessing stable dynamics and not merely as a searching 
for stable sub-systems out of which more complicated systems can be built 
up. The fact that more complex living organisms have bigger amounts of 
genetic material, together with the fact that real genetic networks actu- 
ally exhibit high homeostatic stability, suggest that evolutionary processes 
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consist of both kinds of "searching," hierarchical-modular complexification 
and dynamical stability. 

So far we have considered processes occurring at the level of genomes, 
but evolution also acts at the level of populations, making the organisms of 
a given species become better adapted to their environment. Living beings 
are subjected to all kinds of external fluctuations, and the survival of the 
species depends on the capability of its members to recover from those ran- 
dom perturbations. Evolutionary processes produce organisms with a high 
degree of homeostatic stability and of adaptability to the environment, 
even in the presence of external fluctuations. When thinking of evolution- 
ary processes, one usually supposes that living organisms tend to exclude 
noise since stability is more conductive to functioning than chaos, and also 
that evolutionary processes are able to recognize and favor such stability. 
Nevertheless, Michael D. Stern has pointed out that none of these assump- 



tions has been rigorously proven yet ( Stern [1999 1 ) . Furthermore, he has 
shown that, under certain circumstances, noise not only is not excluded 
from an adaptative system, but it is required for the adaptation of the 
system to the environment; without noise, such adaptation would not be 
possible. 

In this work Stern considers a population of M organisms, each of which 
is a N-K model composed of N = 100 elements. Each organism is in the 
critical phase (K = 2 and p = 0.5) and therefore, according to Kauffman, 
they are in the state of highest homeostatic stability. The linkages among 
the elements and the evolution rules fi are assigned in the usual way, but 
additionally, a noisy signal rj(t) is applied to m randomly chosen elements 
of every organism. The noise is introduced through only one of the two 
inputs of each one of the m elements. Suppose for example that and 
<7i 2 are the two inputs of Cj. With no noise, the value of cr, at time t + 1 
would be given, as usual, by 

a i (t+l) = f i (a il (t),a i2 {t)). 

But if <Ji is one of the m elements to which the noisy signal i](t) is being 
applied, then the value of Ui at time t + 1 is given now by 

<j i (t+l) = f i (a il (t),ri(t)). 

Note that the noisy signal rj{t) is the same for the m elements of every 
one of the M organisms. This takes into account the fact that in a real 
population, every organism in the population is subjected to the same 
(noisy) environment. 

The evolution of the population is now determined by an external cri- 
teria which has to be fulfilled. This is the phenotype of the population on 
which natural selection will be acting. In Stern's work, the phenotype to be 
selected is an integer time series, obtained by counting the number of pos- 
itive states (+1) occurring in a given subset of elements (output elements) 
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in each organism. This time series is to fit a predefined time function F(t) 
(the target function) , and in each generation the organisms better adapted 
to F(t) are selected. 

Selection of the organisms is made as follows. At the beginning of each 
generation, every organism is replicated R times allowing some mutations 
which consist mainly of randomly moving one of the input connections of 
a randomly chosen element, and randomly changing the Boolean function 
of a randomly chosen element. So, at the beginning of every generation 
the population actually consists of R x M elements. The system is then 
evolved during 100 time steps, after which the M elements which best fit 
the external criteria F(t) are selected, starting another generation. 
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FIGURE 6.3. (C) Random output of the starting network (solid), the target func- 
tion (dotted), and the fixed binary noise added to the system (lower) in Stern's 
model. (D) The evolved network generates a good approximation of the target 
function after 2500 generations. (E) The same network as in D but operating 
in the absence of noise fails completely in fitting the target function, indicating 
that the evolved network has been "imprinted" by the arbitrary noise sequen ce 
present during its evolution and requires it to function (taken from Stern [1999 ) . 



The results of the simulation are shown in Fig. 6.3, in which three output 
series are compared with the target function F(t) (dotted line). Fig. 3.3G 
shows the output signal of the initial generation, which has not yet passed 
through any selection process. In the presence of noise 77(f), after 2500 gen- 
erations the target function is very well approximated by the output signal 
of the evolving network, as shown if Fig. |6.3| D. However, in the absence of 
noise the same network fails completely to approximate the target function 
(Fig. S.3E). This surprising result implies that noisy perturbations may be 
essential for the adaptation of organisms to their environments. Stern has 
called this phenomenon noise imprinting in evolution, and claims that it 
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"may be a prototype of a general form of symmetry breaking that leads to 
the evolution of structures that are suboptimal and unnecessarily depen- 
dent on fortuitous features of the environment" . 



6.3 Social Systems 

It is always difficult to extrapolate the models and techniques of physics 
to the realm of social sciences. To begin with, it is not yet clear how to 
define a social system in a way that is suitable for the kind of analysis 
that physicists are used to. The complexity we see in human societies is of 
a quite different nature than the one present in physical systems, so one 
must be extremely cautious when talking about "complexity" and "chaos" 
in social sciences. Nonetheless, some parallels can be established between 
the complex behaviors of Boolean networks and human organizations. 

Boolean networks were first introduced in social sciences by Schclling 
in 1971 flShclling [197l| ). Following [Kluver and Schmidt [199^ ], the main 



aspect in a social system is that social actors act according to certain social 
rules; the intended or unintended consequences of their actions generate 
other actions by the same or other actors and so forth. The result of such 
interactions is reflected in the dynamics of the system. 

A key difficulty in the modeling of social networks is determining the 
quantities that convey enough information to describe a given system, and 
which have some degree of predictability. One must also identify the key 
parameters governing the behavior. Two kind of parameters can be distin- 
guished: those which determine the topology of the system, namely, how 
individuals interact, and the ones which determine the social rules, i.e., the 
consequences of the interactions. To illustrate these two kind of parame- 
ters, we can focus our attention on the connectivity K of Boolean networks 
and the weight parameter p associated to the evolution functions fi . 

In ancient societies, characterized by very rigid dictatorial or monarchic 
regimes, the connectivity among people was very small since the majority 
of individuals were restricted to interact only with the few people in their 
communities. Mobility of individuals among different social classes was also 
minimal. Poor people were condemned to remain poor while rich people 
usually remained rich; social mobility, if any, was allow only from upper to 
lower classes. The course of the whole society depended largely on very few 
people (the King and his clique). In contrast, modern democratic societies 
are characterized by much bigger connectivities since individuals have a 
broader spectrum of interactions. In addition, social mobility is greater 
in democratic societies since, at least in principle, anybody can become 
president or actor or any other thing. 

A binary network can be used for the modeling of a society of interacting 
individuals, in which poverty is represented by 1 while richness is repre- 
sented by 0. Of course, in real societies the connectivity among individuals 
varies from one individual to another, and so does the range of influence of 
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their decisions; but it is clear that even if we assign the same connectivity 
to each member of the society, the out-degree will not be the same for all 
individuals^ and therefore the range of influence that each member of the 
society has on the other members varies among them. 

With this setup, a rigid feudal society would be characterized by low 
values of the parameter K and by evolution functions ft whose weight pa- 
rameter p is very low, reflecting the fact that the majority of the individuals 
are poor no matter what they do or whom they interact with. This kind 
of society presents a very simple dynamics in a frozen regime; no changes 
are expected nor unpredictable behavior. The dynamics is governed by a 
very reduced fraction of relevant elements, while the rest of the individuals 
remain in a frozen state with nearly zero mobility. 

On the other hand, a perfect democratic society would be rather charac- 
terized by high connectivities and weighted functions fi whose parameter 
p is close to 0.5. In this situations, we have seen that the system exhibits 
chaotic behavior, with an apparently random dynamics. The society would 
never reach a stable state or a stable attractor, for in this regime the length 
of the cycles is extremely large. 

If these results reflect some of the fundamental aspects behind the dy- 
namics of social organizations, the conclusion would be that the complex- 
ity that we see in modern democratic societies is inherent to the demo- 
cratic principles (parameter values) on which these societies have been 
constructed. Apparently, democracy and complexity are tied together, and 
it is a political decision whether or not it is worthwhile to sacrifice democ- 
racy in the interest of predictability. We should stress, though, that Boolean 
networks are far from accurate representations of a human society, so we 
should consider the previous results only as indicative trends and not as 
the matter of the facts. 



6.4 Neural Networks 

The first neural network model was introduced nearly 60 years ago by 
McCulloch and Pitts i n an attempt to understand some of the cognitive 
processes of the brain ([McCulloch and Pitts [194^] ). The topic of neural 



networks has since grown enormously, and presently it covers a great va- 
riety of fields, ranging from neurophysiology to computational algo rithms 
for pattern recognition and non-parametric optimization processes flChcng 



and Tittcrington [1994 ]) . Nonetheless, it is somehow accepted that the two 
major applications of neural networks consist in the understanding of real 
neural sys tems (the brai n) , and in the development of machine- learning al- 



gorithms ( Farmer [199C| ] ) . In this section, we will address the subject from 



12 The out-degree of node <r; is the number of different nodes which are affected by &i, 
and can run from up to TV. 
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a point of view based on statistical physics. 

As in the N-K model, a typical neural network consists of a set of N 
nodes, each of which receives inputs from K other nodes. The difference 
from the N-K model is that in a neural network the dynamics are given 



cri{t + 1) = Sign V" w il a i] (t) + h 



K 



(6.5) 



where the synapse weights {lUy} and the activation threshold h are random 
variables that can adapt (in learning processes) or stay fixed (in recall or 
association processes). Note that with the above definitions, the nodes have 
the values ai = ±1. 

K. E. Kurten has considered the case in which the synapse weights are 
independent random variables distributed according to a symmetric prob- 
ability function P(ui). By analyzing the Hamming distance, he has shown 
that this kind of neural network has behavior similar to that of the N- 
K mode l, in that the t hree phases (frozen, critical and chaotic) are also 
present ( Kurten [1988a jj]), depending on the connectivity K and the dilu- 



tion of the network (dilution is a measure of the amount of nodes for which 
the synapse weight is 0). The main result is that for low connectivity and 
dilution, the neural network and the N-K model exhibit exactly the same 
dynamics in the Hamming distance. 

Another important application of neural networks consists in investigat- 
ing the tolerance to the influence of external noise in the dynamical orga- 
nization of the network. Huepe and Aldana-Gonzalcz [2002] ] have analyzed 



3 There are other choices for the dynamics, such as 



( K \ 

tTi(t + 1) = tanh ^ WijCTij (t) + h J , 



which produces an output between —1 and +1, or 



<Xi(i + 1) = j 1 + exp ^ ^2 w i ju ij (t) + h 
3=1 



with an output between and 1. Possible choices abound, depending on the particular 
task for which the network has been designed. We will f ocus on the particular class of 



neural networks obeying the dynamics given by equation 6.5 
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FIGURE 6.4. Bifurcation diagram of the order par ame ter as a function of the 
noise intensity rj, for the neural network model in (6.6). In the case shown, the 
network was composed of N = 10 J elements with connectivity K = 11, and the 
synapse weights Wij were uniformly distributed between and 1. The inset shows 
the finite-siz e effects in the numerical results ne ar the critical noise j) c ~ 0.2838. 
(Taken from Huepe and Aldana-Gonzalez [2002].) 



a neural network whose dynamics is given by 

{Sign fejLi w ijVij (*)) with probability 1 - r), 
-Sign fejLi w ij a ij{^)) with probability ??. 

(6.6) 

The dynamics can thus be changed from purely deterministic to purely 
random by varying the noise intensity 77 between and 1/2. In the case 
rj = the network will typically converge to an ordered state in which 
all the oi acquire the same value. For other values of 77, the instantaneous 
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amount of order in the networks is determined by the parameter 



i 



N 



lim — > a; (t) , 

1=1 



(6.7) 



so that \s(t)\ m 1 for an "ordered" system in which most elements take the 
same value, while \s(t)\ « for a "disordered" system in which the elements 
randomly take values +1 or — 1. The order parameter of the network is 
defined then as the time-independent quantity 



* = lim — 



\s(t)\dt. 



(6.8) 



Huepe and Aldana have shown tha t, under very general conditions, the 
neural network model given in (6.6) undergo es a dynamical second order 
phase transition, as it is illustrated in Fig. 3.4 More specifically, if the 
probability density ~P(w) of the synapse weights is a non-symmetric, but 
otherwise arbitrary function, and if the linkages of the network are chosen 
randomly, then there exists a critical value rj c of the noise in the vicinity of 
which the order parameter behaves as 



for r] > Tfc, 



(6.9) 







where C and r\ c are constants whose values depend on the connectivity 
K and on the particular form of P(w). It is worth mentioning that the 
dynamics of the network are exactly the same if both the linkages and 
the synapse weights Wij are time-independent, or if they are re-assigned at 
every time step. As far as we know, this is the first network model in which 
it has been shown that the quenched and annealed approaches coincide 
exactly. 

Other models of particular interest are those in which the synapse weights 
Wij are in turn given by "sub-layers" of connections, such as in the Little- 



Hopfield model (Hopfield [1982] 



Little 



1974]), 



where = ±1 is the value of site i in the sub-layer /i, and the coefficients 
Cij are independent random variables distributed according to some given 
probability distribution P(C). These networks are important because by 
suitable changing the coefficients C»j , the network can be cither synchro- 



nized or "taught" some particular task 


Hopfield [19991; Abarbanel, Rabi-| 


novich, Selverston, and Bazhcnov [1996 


; Dcrrida, Gardner, and Zippcliue 


[1987 


I)- 
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